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Abstract 



We present direct methods, algorithms, and symbohc software for the computation of con- 
servation laws of nonlinear partial differential equations (PDEs) and differential-difference 

equations (DDEs). 

Our method for PDEs is based on calculus, linear algebra, and variational calculus. First, 
we compute the dilation symmetries of the given nonlinear system. Next, we build a candidate 
density as a linear combination with undetermined coefficients of terms that are scaling 
invariant. The variational derivative (Euler operator) is used to derive a linear system for 
the undetermined coefficients. This system is then analyzed and solved. Finally, we compute 
the flux with the homotopy operator. 

The method is applied to nonlinear PDEs in (1 + 1) dimensions with polynomial nonlinear- 
ities which include the Korteweg-de Vries (KdV), Boussinesq, and Drinfel'd-Sokolov- Wilson 
equations. An adaptation of the method is applied to PDEs with transcendental nonlineari- 
ties. Examples include the sine-Gordon, sinh-Gordon, and Liouville equations. For equations 
in laboratory coordinates, the coefficients of the candidate density are undetermined functions 
which must satisfy a mixed linear system of algebraic and ordinary differential equations. 

For the computation of conservation laws of nonlinear DDEs we use a splitting of the 
identity operator. This method is more efficient that an approach based on the discrete 
Euler and homotopy operators. We apply the method of undetermined coefficients to the 
Kac-van Moerbckc, Toda, and Ablowitz-Ladik lattices. To overcome the shortcomings of 
the undetermined coefficient technique, we designed a new method that first calculates the 
leading order term and then the required terms of lower order. That method, which is no 
longer restricted to polynomial conservation laws, is applied to discretizations of the KdV and 
modified KdV equations, and a combination thereof. Additional examples include lattices 
due to Bogoyavlenskii, Belov-Chaltikian, and Blaszak-Marciniak. 

The undetermined coefficient methods for PDEs and DDEs have been implemented in 
Mathematica. The code TransPDEDensityFlux.m computes densities and ffuxes of systems 
of PDEs with or without transcendental nonlinearities. The code DDEDensityFlux.m does 
the same for polynomial nonlinear DDEs. Starting from the leading order terms, the new 
Maple library discrete computes densities and fluxes of nonlinear DDEs. 

The software can be used to answer integrability questions and to gain insight in the phys- 
ical and mathematical properties of nonlinear models. When applied to nonlinear systems 
with parameters, the software computes the conditions on the parameters for conservation 
laws to exist. The existence of a hierarchy of conservation laws is a predictor for complete 
integrability of the system and its solvability with the Inverse Scattering Transform. 
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1 Introduction 

This chapter focuses on symbolic methods to compute polynomial conservation laws of par- 
tial differential equations (PDEs) in (1 + 1) dimensions and differential-difference equations 
(DDEs), which arc semi-discrctc lattices. For the latter we treat systems where time is 
continuous and the spatial variable has been discretized. 

Nonlinear PDEs that admit conservation laws arise in many disciplines of the applied 
sciences including physical chemistry, fluid mechanics, particle and quantum physics, plasma 
physics, elasticity, gas dynamics, electromagnetism, magncto-hydro-dynamics, nonlinear op- 
tics, and the bio-sciences. Conservation laws are fundamental laws of physics that maintain 
that a certain quantity will not change in time during physical processes. Familiar conser- 
vation laws include conservation of momentum, mass (matter), electric charge, or energy. 
The continuity equation of electromagnetic theory is an example of a conservation law which 
relates charge to current. In fluid dynamics, the continuity equation expresses conservation 
of mass, and in quantum mechanics the conservation of probability of the density and flux 
functions also yields a continuity equation. 

There are many reasons to compute conserved densities and fluxes of PDEs explicitly. 
Invariants often lead to new discoveries as was the case in soliton theory. One may want to 
verify if conserved quantities of physical importance (e.g. momentum, energy, Hamiltonians, 
entropy, density, charge) are intact after constitutive relations have been added to close a 
system. For PDEs with arbitrary parameters one may wish to compute conditions on the 
parameters so that the model admits conserved quantities. Conserved densities also facilitate 
the study of quahtative properties of PDEs [86, 97], such as recursion operators, bi- or 
tri-Hamiltonian structures, and the like. They often guide the choice of solution methods or 
reveal the nature of special solutions. For example, an infinite sequence of conserved densities 
is a predictor of the existence of solitons [7] and complete integrability [2] which means that 
the PDE can be solved with the Inverse Scattering Transform (1ST) method [2] . 

Conserved densities aid in the design of numerical solvers for PDEs [87, 88] and their sta- 
bility analysis (see references in [23]). Indeed, semi-discretizations that conserve discrete con- 
served quantities lead to stable numerical schemes that are free of nonlinear instabilities and 
blowup. While solving DDEs, which arise in nonlinear networks and as semi-discretizations 
of PDEs, one should check that their conserved quantities indeed remain unchanged as time 
steps are taken. 
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Computer algebra systems (CAS) like Mathematica, Maple, and REDUCE, can greatly assist 
the computation of conservation laws of nonlinear PDEs and DDEs. Using CAS interactively, 
one can make a judicious guess (ansatz) and find a few simple densities and fluxes. Yet, 
that approach is fruitless for complicated systems with nontrivial conservation laws with 
increasing complexity. Furthermore, completely integrable equations PDEs [2, 7, 74, 89] 
and DDEs [10, 75] admit infinitely many independent conservation laws. Computing them 
is a challenging task. It involves tedious computations which are prone to error if done 
with pen and paper. Kruskal and collaborators demonstrated the complexities of calculating 
conservation laws in their seminal papers [67, 78, 79] on the Korteweg-de Vries (KdV) equation 
from sohton theory [2, 7, 30]. We use this historical example to introduce the method of 
undetermined coefficients. 

In the first part of this chapter we cover the symbolic computation of conservation laws 
of completely integrable PDEs in (1 + 1) dimensions (with independent variables x and t). 
Our approach [8, 38, 51, 53] uses the concept of dilation (scahng) invariance and the method 
of undetermined coefficients. Our method proceeds as follows. First, build a candidate den- 
sity as a linear combination (with undetermined coefficients) of "building blocks" that are 
homogeneous under the scaling symmetry of the PDE. If no such symmetry exists, construct 
one by introducing parameters with scaling. Next, use the Euler operator (variational deriva- 
tive) to derive a linear algebraic system for the undetermined coefficients. After the system 
is analyzed and solved, use the homotopy operator to compute the flux. When applied to 
systems with parameters, our codes can determine the conditions on the parameters so that 
a sequence of conserved densities exists. 

The method is applied to nonlinear PDEs in (1-1-1) dimensions with polynomial terms 
which include the KdV, Boussinesq, and Drinfel'd-Sokolov- Wilson equations. An adaptation 
of the method is applied to PDEs with transcendental nonlinearities. Examples include 
the sine-Gordon, sinh-Gordon, and Liouville equations. For equations written in laboratory 
coordinates, the coefficients of the candidate density are undetermined functions which must 
satisfy a mixed linear system of algebraic and ordinary differential equations (ODEs). 

Capitalizing on the analogy between PDEs and DDEs, the second part of this chapter deals 
with the symbohc computation of conservation laws of nonlinear DDEs [34, 39, 42, 51, 54, 57]. 
Again, we use scaling symmetries and the method of undetermined coefficients. One could 
use discrete versions of the Euler operator (to verify exactness) and the homotopy operator 
(to invert the forward difference). Although these operators might be valuable in theory, they 
are highly inefficient as tools for the symbolic computation of conservation laws of DDEs. We 
advocate the use of a "splitting and shifting" technique, which allows us to compute densities 
and fluxes simultaneously at minimal cost. The undetermined coefficient method for DDEs 
is illustrated with the Kac-van Moerbeke, Toda, and Ablowitz-Ladik lattices. 

There is a fundamental difference between the continuous and discrete cases in the way 
densities are constructed. The total derivative has a weight whereas the shift operator does 
not. Consequently, a density of a PDE is bounded in order with respect to the space variable. 
Unfortunately, there is no a priori bound on the number of shifts in the density, unless a 
leading order analysis is carried out. To overcome this difficulty and other shortcomings 
of the undetermined coefficient method, we present a new method to compute conserved 
densities of DDEs. That method no longer uses dilation invariance and is no longer restricted 
to polynomial conservation laws. Instead of building a candidate density with undetermined 
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coefficients, one first computes the leading order term in the density and, secondly, generates 
the required terms of lower order. The method is fast and efficient since unnecessary terms 
are never computed. The new method is illustrated using a modified Volterra lattice as 
an example. The new method performs exceedingly well when applied to lattices due to 
Bogoyavlenskii, Belov-Chaltikian, and Blaszak-Marciniak. The new method is also apphed 
to completely integrable discretizations of the KdV and modified KdV (mKdV) equations, 
and a combination thereof, known as the Gardner equation. Starting from a discretized 
eigenvalue problem, we first derive the Gardner lattice and then compute conservation laws. 

There are several methods (see [51]) to compute conservation laws of nonhnear PDEs 
and DDEs. Some methods use a generating function [2, 7], which requires the knowledge of 
key pieces of the 1ST. Another common approach uses the link between conservation laws 
and symmetries as stated in Noether's theorem [14, 15, 66, 81]. However, the computation 
of generalized (variational) symmetries, though algorithmic, is as daunting a task as the 
direct computation of conservation laws. Most of the more algorithmic methods [12, 13, 20, 
25, 63, 101], require the solution of a determining system of ODEs or PDEs. Despite their 
power, only a few of these methods have been implemented in CAS. We devote a section to 
symbolic software for the computation of conservation laws. Additional reviews can be found 
in [38, 51, 101]. 

Over the past decade, in collaboration with students and researchers, we have designed and 
implemented direct algorithms for the computation of conservation laws of nonlinear PDEs 
and DDEs. We purposely avoid Noether's theorem, pre-knowledge of symmetries, and a La- 
grangian formulation. Neither do we use differential forms or advanced differential-geometric 
tools. Instead, we concentrate on the undetermined coefficient method for PDEs and DDEs, 
which uses tools from calculus, hnear algebra, and the variational calculus. Therefore, the 
method is easy to implement in Mathematica and easy to use by scientists and engineers. 
The code TransPDEDensityFlux.m computes densities and fluxes of systems of PDEs with 
or without transcendental nonlinearities. The code DDEDensityFlux.m does the same for 
polynomial nonlinear DDEs. Starting from the leading order terms, the new Maple library 
discrete computes densities and fluxes of nonlinear DDEs very efficiently. The software 
can thus be used to answer integrability questions and to gain insight in the physical and 
mathematical properties of nonlinear models. 

Our software is in the public domain. The Mathematica packages and notebooks are 
available at [48] and Hickman's code in Maple is available at [56]. We are currently working on 
a comprehensive package to compute conservation laws of PDEs in multiple space dimensions 
[45, 51, 83]. 

Part I: Partial Differential Equations in (1+1) Dimensions 

In this flrst part we cover PDEs in (1 + 1) dimensions, that is, PDEs in one space variable 
and time. Starting from a historical example, we introduce the concept of dilation invariance 
and use the method of undetermined coefficients to compute conservation laws of evolution 
equations. Later on, we adapt the method of undetermined coefficients to cover PDEs with 
transcendental terms. 
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2 The Most Famous Example in Historical Perspective 



The story of conservation laws for nonlinear PDEs begins with the discovery of an infinite 
number of conservation laws of the ubiquitous Korteweg-de Vries equation which models a 
variety of nonlinear wave phenomena, including shallow water waves [46] and ion-acoustic 
waves in plasmas [2, 7, 30]. The KdV equation can be recast in dimensionless variables as 



Ut + aUUx + U2,a: = 0, 



where the subscripts denote partial derivatives, i.e. Ut = ^,Ux = and u^x = f^- The 
parameter a can be scaled to any real number. Commonly used values are a = ±1 or a = ±6. 
Equation (1) is an example of a scalar (1 + 1)— dimensional evolution equation, 

Ut^ F{x,t,U,U^,U2x,- ■ ■ ,U'nx), (2) 

of order n in the independent space variable x and of first order in time t. Obviously, the 
dependent variable is u{x, t). If parameters are present in (2), they will be denoted by lower- 
case Greek letters. A conservation law of (2) is of the form 

Dip+D^J = 0, (3) 

which is satisfied for all solutions u{x,t) of the PDE. In physics, p is called the conserved 
density (or charge); J is the associated flux (or current). In general, both are differential 
functions (functionals) , i.e. functions of x,t,u, and partial derivatives of u with respect to x. 
In (3), D-r denotes the total derivative with respect to x, that is, 

^ ^ dJ dJ 

fe=0 

where N is the order of J, and is the total derivative with respect to t, defined by 

o„.|.,>,H|.f;^oj.. (5, 

k=0 

where p'[ut] is the Frechet derivative of p in the direction of Ut and M is the order of p. 

The densities p^-^^ = u and p^^^ = of (1) were long known. In 1965, Whitham [98] had 
found a third density, p^^^ = — ^u^, which, in the context of water waves, corresponds to 
Boussinesq's moment of instability [76] . One can readily verify that 

Dt{u) + D,{^au' + U2x) = 0, (6) 

Dt{u'') + Dx{^au^ -ul + 2uu2x)^0, (7) 

3 3 3 6 

Dt{u^ ul) + D^{-au^ - 6uul + 3u\2x + -«L ^^x^sx) = 0. (8) 

q; 4 a a 
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Indeed, (6) is the KdV equation written as a conservation law; (7) is obtained after multiplying 
(1) by 2u; (8) requires more work. Hence, the first three density-fiux pairs of (1) are 

^au'^ + U2x, (9) 
2 

-au^ -ul + 2uu2x, (10) 

^au^ - Guul + 3u'^U2x + — ~ —UxUsx- (H) 
4 a a 

Integrals of motion readily follow from the densities. Indeed, assuming that J vanishes at 
infinity (for example due to sufficiently fast decay of u and its x derivatives), upon integration 
of (3) with respect to x one obtains that 

/oo 
pdx (12) 
-oo 

is constant in time. Such constants of motion also arise when u is periodic, in which case one 
integrates over the finite period. Depending on the physical setting, the first few constants 
of motion (i.e. integrals (12)) express conservation of mass, momentum, and energy. 

Martin Kruskal and postdoctoral fellow Norman Zabusky discovered the fourth and fifth 
densities for the KdV equation [111]. However, they failed in finding a sixth conservation 
law due to an algebraic mistake in their computations. Kruskal asked Robert Miura, also 
postdoctoral fellow at the Princeton Plasma Physics Laboratory at New Jersey, to search for 
further conservation laws of the KdV equation. Miura [78] computed the seventh conservation 
law. After correcting the mistake mentioned before, he also found the sixth and eventually 
three additional conservation laws. Rumor [80] has it that in the summer of 1966 Miura went 
up into the Canadian Rockies and returned from the mountains with the first 10 conservation 
laws of the KdV equation engraved in his notebook. This biblical metaphor probably does 
not do justice to Miura's intense and tedious work with pen and paper. 

With ten conservation laws in hand, it was conjectured that the KdV equation had an 
infinite sequence of conservation laws, later proven to be true [67, 79]. Aficionados of explicit 
formulas can find the first ten densities (and seven of the associated fiuxes) in [79] and 
the eleventh density (with 45 terms) in [67], where a recursion formula is given to generate 
all further conserved densities. As an aside, in 1966 the first five conserved densities were 
computed on an IBM 7094 computer with FORMAC, an early CAS. The sixth density could 
no longer be computed because the available storage space was exceeded. In contrast, using 
a method of undetermined coefficients, the first eleven densities were computed in 1969 on 
a AEC CDC-6600 computer in a record time of 2.2 seconds. Due to hmitations in handhng 
large integers, the computer could not correctly produce any further densities. 

Undoubtedly, the discovery of conservation laws played a pivotal role in the comprehensive 
study of the properties and solutions of nonlinear completely integrable PDEs (like the KdV 
equation) and the development of the 1ST (see e.g. [80] for the history). Clifford Gardner, 
John Greene, Martin Kruskal, and Robert Miura received the 2006 Leroy P. Steele Prize [115], 
awarded by the American Mathematical Society, for their seminal contribution to research 
on the KdV equation. In turn, Martin Kruskal has received numerous honors and awards 
[114] for his fundamental contributions to the understanding of integrable systems and soliton 
theory. This chapter is dedicated to Martin Kruskal (1925-2006). 



p(3) ^ ^3 _ 3^2^ ^(3) _ 
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3 The Method of Undetermined Coefficients 



We now sketch the method of undetermined coefficients to compute conservation laws [38, 96], 
which draws on ideas and observations in before mentioned work by Kruskal and collaborators. 



3.1 Dilation Invariance of Nonlinear PDEs 

Crucial to the computation of conservation laws is that (3) must hold on the PDE. This is 
achieved by substituting ut (and Utx,Utxx, etc.) from (1) in the evaluation of (6)-(8) and in 
all subsequent conservation laws of degree larger than 3. The elimination of all derivatives 
of u in favor of x derivatives has two important consequences: (i) any symmetry of the PDE, 
in particular, the dilation symmetry, will be adopted by the conservation law, (ii) once Dj is 
computed and evaluated on the PDE, t becomes a parameter in the computation of the flux. 

We will first investigate the dilation (scaling) symmetry of evolution equations. The KdV 
equation is dilation invariant under the scaling symmetry 

{t, X, u) {X-H, X-^x, X\) , (13) 

where A is an arbitrary parameter. Indeed, after a change of variables with i — X''H, x — 
X~^x,u = X'^u, and cancellation of a common factor A^, the KdV for u{x,i) arises. The 
dilation symmetry of (1) can be expressed as 

dx'^ ' dt dx^ ' 

which means that u corresponds to two x— derivatives and the time derivative corresponds to 
three derivatives. If we define the weight, W, of a variable (or operator) as the exponent 
of A in (13), then W{x) = -1 or W{-§;^) = 1; W{t) = -3 or W{-§-^) = 3, and W{u) = 2. 

All weights of dependent variables and the weights of d/dx,d/dt, are assumed to be 
non-negative and rational. The rank of a monomial is defined as the total weight of the 
monomial. Such monomials may involve the independent and dependent variables and the 
operators D^-, ^, and D^. Ranks must be positive integers or positive rational numbers. 

An expression (or equation) is uniform in rank if its monomial terms have equal rank. 
For example, (1) is uniform in rank since each of the three terms has rank 5. 

Conversely, if one does not know the dilation symmetry of (1), then it can be readily 
computed by requiring that (1) is uniform in rank. Indeed, setting W{d/dx) — 1 and equating 
the ranks of the three terms in (1) gives 

d 

W{u) + W{—) = 2W{u) + l = W{u) + 3, (15) 

which yields W{u) — 2,W{d/dt) — 3, and, in turn, confirms (13). So, requiring uniformity 
in rank of a PDE allows one to compute the weights of the variables (and thus the scaling 
symmetry) with linear algebra. 

Dilation symmetries, which are special Lie-point symmetries, are common to many non- 
linear PDEs. Needless to say, not every PDE is dilation invariant, but non- uniform PDEs can 
be made uniform by extending the set of dependent variables with auxiliary parameters with 
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appropriate weights. Upon completion of the computations one can set these parameters to 
one. In what follows, we set W{d/dx) = 'I^(D^) = 1 and W{d/dt) = W{Dt). Applied to 
(6), rankp*^^'' = 2, rankJ'-^-* = 4. Hence, rank (Df p*^^'') = rank (D-j, J*^^'') = 5. Therefore, (6) 
is uniform of rank 5. In (7), rankp^^^ = 4 and rank J^^^ = 6, consequently, (7) is uniform of 
rank 7. 

In (11), each term in p*^'^^ has rank 6 and each term in J^^^ has rank 8. Consequently, 
rank (D^ p^'^^) = rank (D^; J'-'^-') = 9, which makes (8) is uniform of rank 9. All densities of (1) 
are uniform in rank and so are the associated fluxes and the conservation laws. 

Equation (1) also has density-flux pairs that depend explicitly on t and x; for example, 

2-2 / 2 \ 2 

p — tu^ H — xu, J — t{-au^ — + 2uu2x) — x [u^ U2x ) H — u^- (16) 

a 3 \ a J a 

Since W{x) = —1 and W{t) = —3, one has rankp = 1, and rank J = 3. The methods 
and algorithms discussed in subsequent sections have been adapted to compute densities and 
fluxes explicitly dependent on x and t. Instead of addressing this issue in this chapter, we 
refer the reader to [38, 53]. 

3.2 The Method of Undetermined Coefficients AppHed to a Scalar 
Nonhnear PDE 

We outline how densities and fluxes can be constructed for a scalar evolution equation (2). 
To keep matters transparent, we illustrate the steps for the KdV equation resulting in p^^^ of 
rank R — 6 with associated flux J^^^ of rank 8, both hsted in (11). The tools needed for the 
computations will be presented in the next section. 

• Select the rank R of p. Make a list, TZ, of all monomials in u and its x-derivatives so 
that each monomial has rank R. This can be done as follows. Starting from the set V of 
dependent variables (including parameters with weight, when applicable), make a set Ai of 
all non-constant monomials of rank R or less (but without x'— derivatives). Next, for each 
term in A4, introduce the right number of x— derivatives to adjust the rank of that term. 
Distribute the x— derivatives, strip off the numerical coefficients, and gather the resulting 
terms in a set TZ. For the KdV equation and R = 6, V = {u} and Ai — {u^,u^,u}. Since 
u^jU^, and u have ranks 6,4 and 2, respectively, one computes 

Ignoring numerical coefficients in the right hand sides of the equations in (17), one gets 

TZ = {u^,ul,UU2x,U4x}- 

• Remove from TZ all monomials that are total x— derivatives. Also remove all "equivalent" 
monomials, i.e. the monomials that differ from another by a total x— derivative, keeping the 
monomial of lowest order. Call the resulting set S. In our example, u^^x must be removed 
(because U4x = DxUsx) and uu2x must be removed since uu2x and are equivalent. Indeed, 
uu2x — ^x{uux) — ul- Thus, S = {u^,ul}. 

• Linearly combine the monomials in S with constant undetermined coefficients Cj to obtain 
the candidate p. Continuing with the example, 

p^ciu^ + C2ul, (18) 
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which is of first order in x. 

• Using (5), compute p. Apphed to (18) where M — 1, one gets 

Dtp^ (3ciii'l + 2c2U^r>^)[ut]. (19) 
As usual, D° = I is the identity operator. 

• Evaluate — D^p on the PDE (2) by replacing Ut by F. The result is a differential function 
E in which i is a parameter. For the KdV equation (1), F = —{auux + Us^)- After reversing 
the sign, the evaluated form of (19) is 

E = {3CiU^\ + 2C2Ua:D^){aUU^ + U3a:) 

— SCiaU^Ux + 2C20!ul + 2C2aUUxU2x + "iCiU^U^x + 2C2UxU4x- (20) 

• To obtain a conservation law, E must be a total derivative. Starting with highest orders, 
repeatedly integrate E by parts. Doing so, allows one to write E as the sum of a total 
x— derivative, D^, J, and a non-integrablc part (i.e. the obstructing terms). J is the (candidate) 
flux with rank J — R-\-W{Dt) — 1. Integration by parts of (20) gives 

3 

E = Dj;(-CiQ;M^ + 3CiU^U2x + C2auul + 2c2UxUsx) - GCiUUxU2x - 2c2U2xUzx + C2au\ 

3 

= Dj,(-ciq;m^ - 3ciuul + C2auul + 3ciu'^U2x + 2c2UxU3x - C2u\^) + (3ci + C2a)ul. (21) 

The candidate flux therefore is 
3 

J — -Ciau'^ — SCiUU^. + C20tUu1. + 3CiU^U2x + 2c2UxU2x — C2ulx- (22) 

• Equate the coefficients of the obstructing terms to zero. Solve the linear system for the 
undetermined coefficients q. In the example (3ci + C2a)u^ is the only obstructing term which 
vanishes for C2 = ~f Ci, where Ci is arbitrary. 

• Substitute the coefficients q into p and J to obtain the final forms of the density and 
associated flux (with a common arbitrary factor which can be set to 1). Setting ci — 1 and 
substituting C2 = — | into (18) and (22) yields p^^^ and J^^^ as given in (11). 

Constructing "minimal" densities, i.e. densities which are free of equivalent terms and total 
derivatives terms, becomes challenging if the rank R is high. Furthermore, integration by 
parts is cumbersome and prone to mistakes if done by hand. Moreover, it would be advan- 
tageous if the integration by parts could be postponed until the undetermined coefficients q 
have been computed and substituted in E. Ideally, the computations of the density and the 
flux should be decoupled. There is a need for computational tools to address these issues, in 
particular, if one wants to compute conservation laws of systems of evolution equations. 

4 Tools from the Calculus of Variations and Differential 
Geometry 

A scalar differential function E of order M is called exact (integrable) if and only if there 
exists a scalar differential function J of order M — 1 such that E — Dx J. Obviously, J = 
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D~^E = J E dx is then the primitive (or integral) of E. Two questions arise: (i) How can 
one test whether or not E is exact? (ii) If E is exact, how can one compute J without 
using standard integration by parts? To answer the first question we wiU use the variational 
derivative (Euler operator) from the calculus of variations. To perform integration by parts 
we will use the homotopy operator from differential geometry. 



4.1 The Continuous Variational Derivative (Euler Operator) 



The continuous variational derivative, also called the Euler operator of order zero, , for 



variable u{x) is defined [81] by 



'(0) _ V"/ \kdE 



fc=0 

dE ^ dE ^2 dE . dE , dE 

du du^ du2a, du^^ duux 

where is a differential function in u{x) of order M. 

A necessary and sufficient condition for a differential function E to be exact is that 
C^^lx)-^ = 0. A proof of this statement is given in e.g. [67]. If Cfl^^-^E ^ 0, then E is not a 
total X— derivative due to obstructing terms. 

Application 1. Returning to (20), we now use the variational derivative to determine Ci 
and C2 so that E of order M — A will be exact. Using nothing but differentiations, we readily 
compute 



«(^) du ""du^ ''du2x ""dus^ ""du^^ 

— 9ciau^Ux + 2c2aUxU2x + Ociw-Ua^ — ^^{'iciau^ + 6c2Q;m^ + 2c2auu2x + 2c2M4a;) 
+ Dl{2c2auu.,) - Dl{?>ciu^) + D^.(2c2Mx) 

= {'dCiaU^Ux + 2C20iUxU2x + QCiUU^x) — {'dCiaU^Ux + l'iC20iUxU2x + 2C2a'UM3x 
+2c2U5x) + {QC2aUxU2x + 2c2aUU2,x) - {l^CiUxU2x + QCiUU^x) + (2C2M5a;) 

= -6(3ci + C2a)uxU2x- (24) 

Note that the terms in v?Ux,uu2,x-, and Ur,^. dropped out. Hence, requiring that £^°^^^£' = 
leads to 3ci+C2Q; = 0. Substituting ci = l,C2 = — |, into (18) yields p^^^ in (11). 

Application 2. It is paramount that the candidate density is free of total derivatives and 
equivalent terms. If such terms were present, they could be moved into the flux J, and their 
coefficients Cj would be arbitrary, p*^"*^-* and p*-^-*, are equivalent if and only if p'^^^ + kp^'^^ = J, 
for some J and non-zero scalar k. We write p^^^ = p^'^\ Clearly p is equivalent to any non-zero 
multiple of itself and p = if and only if p is exact. 

Instead of working with different densities, we investigate the equivalence of terms ti in 
the same density. For example, returning to the set TZ = {u^,ul,uu2x,U4x}, terms ^2 — 
and ts = uu2x are equivalent because ^3 + ^2 = uu2x + = Dx{uUx). 

The variational derivative can be used to detect equivalent and exact terms. Indeed, 
note that vi — £^^^^-^{uu2x) — 2ii2s and V2 — C^^^^^u^ ~ —2u2x are linearly dependent. Also, 
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for ^4 = = ^xU-sx one gets ^3 = C^^^^^s^u^x = 0. To weed out the terms ti in TZ that are 
equivalent or total derivatives, it suffices to check the linear independence of their images Vi 
under the Euler operator. 

One can optimize this procedure by starting from a set TZ where some of the equivalent 
and total derivatives terms have been removed a priori. Indeed, in view of (17), one can 
ignore the highest-order terms (typically the last terms) in each of the right hand sides. 
Therefore, TZ = {u^,u^} and, for this example, no further reduction would be necessary. 
Various algorithms are possible to construct minimal densities. Details are given in [38, 51]. 



4.2 The Continuous Homotopy Operator 

We now discuss the homotopy operator [12, 13, 52, 81] which will allow one to reduce the 
computation of J = D~^E = J E dx to a single integral with respect to an auxiliary variable 
denoted by A (not to be confused with A in Section 3.1). Hence, the homotopy operator 
circumvents integration by parts and reduces the inversion of to a problem of single- 
variable calculus. 

The homotopy operator [81, p. 372] for variable u{x), acting on an exact expression E of 
order M, is given by 

nuix)E = {hE)[Xu]—, (25) 
where the integrand luE is given by 

/„£ = f;fX:M-DJ'-<«')|f. (26) 



k=l \i=0 



kx 



In (25), (/„£') [Alt] means that in one replaces u Xu, — > Aw-r, etc. This is a special 



„,0 „,0 „,0 



2x^ ) "Ma; 



case of the homotopy, A(m*^^) — n'^^^) + between two points, u^^'> — (m°,-u°,-u' 
and u^^^ = {u^.ul, u^^, • • • , wjy^^,), in the jet space. For our purposes we set u^^^ = (0, 0, ■ ■ ■ , 0) 
and u^^^ = {u,Ux,U2x, ■ ■ ■ ,umx)- Formula (26) is equivalent to the one in [52], which in turn 
is equivalent to the formula in terms of higher Euler operators [45, 51]. 

Given an exact differential function E of order M one has J = D^^-E = J E dx — Hu{x)E. 
A proof of this statement can be found in [52]. 

Application. After substituting ci — 1 and C2 — into (20) we obtain the exact expression 

6 

E = Sau^x ~ — 6uUxU2x + 'iu^Usx UxU^x, (27) 

a 

of order M — A. First, using (26), we compute 
4 

luE 




dE 

uDl - UxDx + U2x\)i^ — ) 

OUsx 

dE 

+ {-uDl + UxDl-U2xDx + U3x\){-Q^). (28) 
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After substitution of (27), one gets 



luE = {u\){3au^ + ISul — Quu2x u^x) + {—uD^ + Ux\){—QuUx) 



a 



+ {uDl - Ux\^x + U2x\){'iu^) + {-UDI + Ux\^l - U-2x^x + '"3a:l)( 



6 



Ux) 




a a 



(29) 



which has the correct terms of J'^^^ but incorrect coefficients. Finally, using (25), 



J = nuix)E= / {IuE)[Xu] — 




— -jOiU^ — ^UU^. + Zu^U^x H UxUz: 



/ ( SaA^ii^ — ISA^iiii^ + 9)<^u^U2x H — AmL \uxU2,x I d\ 

V q; q; y 



(30) 



which matches J*^'^^ in (11). 

The crux of the homotopy operator method [12, 13, 26, 81] is that the integration by 
parts of a differential expression like (27), which involves an arbitrary function u{x) and its 
X— derivatives, can be reduced to a standard integration of a polynomial in A. 



5 Conservation Laws of Nonlinear Systems of Polyno- 
mial PDEs 



Thus far we have dealt with the computation of dcnsity-fiux pairs of scalar evolution equa- 
tions, with the KdV equation as the leading example. In this section we show how the 
method and tools can be generalized to cover systems of evolution equations. We will use 
the Drinfel'd-Sokolov- Wilson system and the Boussinesq equation to illustrate the steps. 

5.1 Tools for Systems of Evolution Equations 

For differential functions (like densities and fluxes) of two dependent variables (m, v) and their 
x— derivatives, the total derivatives are 



where Mi and M2 are the (highest) orders of u and v in p, and A''i and A^2 are the (highest) 
orders of u and v in J. 




(31) 




(32) 
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To accommodate two dependent variables, we need Euler operators C:,'s and CiK for 
each dependent variable separately. For brevity, we will use vector notation, that is, i^^^^^-^E — 
(>C^°^^^£',£^°^^£'). Likewise, the homotopy operator in (25) must be replaced by 

n^{x)E = (luE + hE)[Xn] (33) 

where 

Ml /fc-l \ ^ 

/.B^E E-^-eDJ'-'-' £. (34) 

k=l \i=0 J '^^ 

and 

(35) 

k=l \i=0 / '^^ 

where Mi, M2 are the orders of E in li, respectively. In (33), u Xu, Ux — > Xux, ■ ■ ■ , f — > 
Xv, Vx —>■ Xvx, etc. 

5.2 The Drinfel'd-Sokolov- Wilson System: Dilation Invariance and 
Conservation Laws 

We consider a parameterized family of the Drinfel'd-Sokolov- Wilson (DSW) equations 

Ut + SvVx — 0, Vt + 2uVx + aUxV -\- 2vzx — 0, (36) 

where a is a nonzero parameter. The system with a = 1 was first proposed by Drinfel'd and 
Sokolov [31, 32] and Wilson [99]. It can be obtained [59] as a reduction of the Kadomtsev- 
Petviashvili equation (i.e. a two-dimensional version of the KdV equation) and is a completely 
integrablc system. In [109], Yao and Li computed conservation laws of (36), where they had 
introduced four arbitrary coefficients. Using scales on a;, t, u and f , all but one coefficients 
in (36) can be scaled to any real number. Therefore, to cover the entire family of DSW 
equations it suffices to leave one coefficient arbitrary, e.g. a in front of UxV. 

To compute the dilation symmetry of (36), we assign weights, W{u) and W{v), to both 
dependent variables and express that each equation separately must be uniform in rank (i.e. 
the ranks of the equations in (36) may differ from each other). 

For the DSW equations (36), one has 

W{u) + W{d/dt) = 2W{v) + l, 

W{v) + W{d/dt) = W{u) + W{v) + 1 = W{v) + 3, (37) 

which yields W{u) = W{v) = 2, W{d/dt) = 3. The DSW system (36) is thus invariant under 
the scaling symmetry 

(x, t,u,v)^ ( A~^x, X'H, X^u, X^v) , (38) 
where A is an arbitrary scaling parameter. 
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The first tfiree density-flux pairs for the DSW equations (36) are 

P^'^ = J« = (39) 

= V, J^2) ^2{uv + V2a:), if a = 2, (40) 

= (a - ly + -v'', J(3) = 3{auv^ -vl + ^vv^^), (41) 

Both p*^^) and p*^^^ have rank 2; their fluxes have rank 4. The pair (p*-"*^-*, J*-"'^-') exists for any 
a, whereas (p*^^\ J*^^^) only exists if o; = 2. Density p^^"^ of rank 4 and flux J*^^'' of rank 6 are 
valid for any a. At rank R = 6, 

p(^) = (a + l)(a - 2^ - + l)uv' - ^{a - 2)4 - (42) 

The corresponding flux (not shown) has 7 terms. At rank R = 8, 

(5) 4 9 2 2 27 4 9 2 I 3 2 I 45 , oT 2 §1 2 //loN 

provided a ~ 1. The corresponding flux (not shown) has 15 terms. There exists a density-flux 
pair for all even ranks i? < 10 provided a — 1, for which (36) is completely integrable. 

5.3 Computation of a Conservation Law of the Drinfel'd-Sokolov 
Wilson System 

To illustrate how the presence of a parameter, like a, affects the computation of densities, we 
compute p^^^ and p*^^^ of rank R = 2 given in (39) and (40). 

Step 1: Construct the form of the density 

The set of dependent variables is V = {u, v}. Both elements are of rank 2 so, no x— derivatives 
are needed. Thus, A4 = 71 = S = {u,v}. Linearly combining the elements in S gives 

p^ CiU + C2V. 

Step 2: Compute the undetermined coefficients q 

Evaluating E — — D^p = —{ciUt + C2Vt) on (36), yields 

E = Scivv^ + C2{2uv^ + au^v 2v2,x), (44) 

which will be exact if C^^^^^E = {6^^^^E , Cf^^^E) = (0,0). Since E is of order Mi ^ 1 in u 
and order M2 = 3 in v, 

dE dE 

^u{x)^ ^ ~du~ ^''du' ^ ^"^^^^ ~ D^{c2av) = (2 - a)c2V^, (45) 

and 

^fo) ^ dE ^ dE ^2 dE dE 



= ?,ciVx + C2aux — ^x{^civ ^-2c2u) -Dl.{2c2) = {a -2)c2Ua:. (46) 
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Both (45) and (46) will vanish identically if and only if {a — 2)c2 = 0. This equation (with 
unknowns Ci and C2) is parameterized by a 7^ 0. The solution algorithm [38] considers all 
branches of the solution and possible compatibility conditions. Setting ci = 1, leads to either 
(i) C2 = if a 7^ 2, or (ii) C2 arbitrary if a = 2. Setting C2 — 1 leads to the compatibility 
condition, a — 2, and ci arbitrary. Substituting the solutions into p — ciu + C2V gives p — u 
which is valid for any a; and p = -u + C2V or p = CiU + v provided a = 2. In other words, 
p*^^^ = 'U is the only density of rank 2 for arbitrary values of a. For a — 2 there exist two 
independent densities, p^-^^ = u and p^'^^ = v. 

Step 3: Compute the associated flux J 

As an example, we compute the flux in (40) associated with p^"^^ = v and o; = 2. In this case, 
Ci = 0, C2 = 1, for which (44) simplifles into 

E = 2{uvx + u^v + V3x) , (47) 

which is of order Mi — 1 in u and order M2 — 3 in v. Using (34) and (35), we obtain 

luE = u\— = u\{2v) = 2uv, (48) 



and 

OE QE OE 

= {v\){2u) + {vDl - v^D^ + V2j){2) = 2uv + 2v2x. (49) 
Hence, using (33), 

{I^E + hE)[Xvi] y = (4Ak^; + 2^;2^) d\ ^2{uv + V2^), (50) 

which is J*^^-* in (40). The integration of (47) could easily be done by hand. The homotopy 
operator method pays off if the expression to be integrated has a large number of terms. 



5.4 The Boussinesq Equation: Dilation Invariance and Conserva- 
tion Laws 

The wave equation, 

U2t - U2x + + 3'Wti2x + 0CU4x = 0, (51) 

for u{x,t) with real parameter a, was proposed by Boussinesq to describe surface waves in 
shallow water [2]. For what follows, we rewrite (51) as a system of evolution equations, 

ut + Vx- = 0, vt + - 3uux - au3x = 0, (52) 

where v{x,t) is an auxiliary dependent variable. 

The Boussinesq system (52) is not uniform in rank because the terms and auzx lead 
to an inconsistent system of weight equations. To circumvent the problem we introduce an 
auxiliary parameter (3 with (unknown) weight, and replace (52) by 

Ut + Vx^ 0, vt + (5ux - Suux - ausx = 0. (53) 
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Requiring uniformity in rank, we obtain (after some algebra) 

W{u) = 2, W{v) = 3, W{P) = 2, W{^^) = 2. (54) 

Therefore, (53) is invariant under the scahng symmetry 

{x, t, u, V, (5) (A-^x, \-H, X\ Xh, \^(5). (55) 

As the above example shows, a PDE that is not dilation invariant can be made so by extending 
the set of dependent variables with one or more auxiliary parameters with weights. Upon 
completion of the computations one can set each of these parameters equal to 1. 

The Boussinesq equation (51) has infinitely many conservation laws and is completely 
integrable [2, 7]. The first four density-fiux pairs [8] for (53) are 

p(^) = u, J(^) = V, (56) 

p^'^ = V, J(') ^(5u- ^u'' - au2:,, (57) 

^ j(3) ^ ^^^2 _ ^3 _^ 1 ^2 _^ ^aul. - aUU2x, (58) 

= -u^ + y^ + aul, J(^) = 2Puv - 3u% + 2aUxVx - 2au2xV. (59) 

These densities are of ranks 2,3,5 and 6, respectively. The corresponding fiuxes are of one 
rank higher. After setting P — 1 we obtain the conserved quantities of (52) even though 
initially this system was not uniform in rank. 

5.5 Computation of a Conservation Law for the Boussinesq System 

We show the computation of p^^^ and J^^^ of ranks 6 and 7, respectively. The presence of 
the auxiliary parameter f3 with weight complicates matters. At a fixed rank R, conserved 
densities corresponding to lower ranks might appear in the result. These lower-rank densities 
are easy to recognize for they are multiplied with arbitrary coefficients q. Consequently, when 
parameters with weight are introduced, the densities corresponding to distinct ranks are no 
longer linearly independent. As the example below will show, densities must be split into 
independent pieces. 

Step 1: Construct the form of the density 

Augment the set of dependent variables with the parameter j3 (with non-zero weight). Hence, 
V = {u,v,P}. Construct M. = {P'^u, Pu"^, Pu, Pv,u^,u'^,u,v'^,v,uv}, which contains all non- 
constant monomials of (chosen) rank 6 or less (without derivatives). Next, for each term in 
Ai, introduce the right number of ^-derivatives so that each term has rank 6. For example, 

d'^Pu _ d'^u'^ ^2 ^ d^'^ diuv) 

-Q^^Pu2x, -^ = 2u^ + 2uu2x, -^=U4x, —^^^vUx + uVx, etc. (60) 
Gather the terms in the right hand sides of the equations in (60) to get 

7^ = {P^U, Pu^, U^, v'^, VUx, ul, PVx, UVx, pU2x, UU2x, V3x, UAx}- (61) 
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Using (23) and a similar formula for v, for every term ti in TZ we compute Vj = ^u{x)^i ~ 

(^u{x)^'i'^ ^f(x)^i)- ^ (0' 0) ^^^^ discarded and so is Vj. If Vj ^ (0, 0) we verify whether 
or not Vj is linearly independent of the non-zero vectors v^, j — 1, 2, • • • , i — 1. If independent, 
the term ti is kept, otherwise, U is discarded and so is v^. 

Upon application of ^^j,), the first six terms in TZ lead to linearly independent vectors 
Vi through Vfi. Therefore, ti through are kept (and so are the corresponding vectors). For 
tr — (3vx we compute '^7 — ^u{x)('^'^x) — (0, 0). So, tr is discarded and so is V7. For ts — uv^ 

we get V8 = >C^°^^^(iiVa;) = (^;j;, —u^) — — V5. So, is discarded and so is Vg. 

Proceeding in a similar fashion, ^9,^10, in and ti2 are discarded. Thus, TZ is replaced by 

S = {(3\,(3u'^,u^,v^,vu^,ul}, (62) 

which is free of divergences and divergence-equivalent terms. Ignoring the highest-order terms 
(typically the last terms) in each of the right hand sides of the equations in (60) optimizes 
the procedure. Indeed, TZ would have had six instead of twelve terms. Coincidentally, in this 
example no further eliminations would be needed to obtain S. Next, linearly combine the 
terms in S to get 

p = Ci(5'^u -\- C2/9m^ + c^v!^ + c^v^ -\- c^vUx + CqU^. (63) 
Step 2: Compute the undetermined coefficients q 

Compute Dtp. Here, p is of order Mi — 1 in u and order M2 = in v. Hence, application of 
(31) gives 

^ dp, dp ^ dp, 

DtP = —\ut + T;—D^Ut + ^\vt 

OU OUx ov 

= (ci/5^ + 2c2/3m + 3c3U^)ut + {C5V + 2c^u^)utx + (2C4V + Cr,u^)vt. (64) 

Use (53) to eliminate Ut,Utx, and Vt- Then, E — —Dtp evaluates to 

E = (ci/3^ + 2c2(5u 2>czu'^)vx + {c^v 2cqUx)v2x 

^{2ciV C5Ux){l3ux - 3uux - au3x), (65) 

which must be exact. Thus, require that ^^^^^^E = (C^^^^^^E, C^^^^^E) = (0,0). Group like 
terms. Set their coefficients equal to zero to obtain the parameterized system 

/3{c2 - C4) = 0, C3 + C4 = 0, C5 = 0, acs = 0, /3c5 = 0, ac^ - Ce = 0, (66) 

where a ^ and P ^ 0- Investigate the eliminant of the system. Set ci — 1 and obtain the 
solution 

ci = 1, C2 = C4, C3 = -C4, C5 = 0, Cq = ac4, (67) 
which holds without condition on a and f3. Substitute (67) into (63) to get 

p = (3^ + C4 {(3u'^ -u^ + v'^ + aul) . (68) 

The density must be split into independent pieces. Indeed, since C4 is arbitrary, set C4 = or 
C4 = 1, thus splitting (68) into two independent densities, p — (3^u = u and 
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which are p^^^ and p^^^ in (56)- (59). 
Step 3: Compute the flux J 

Compute the flux corresponding to p in (69). Substitute (67) into (65). Take the terms in C4 
and set C4 — 1. Thus, 

E — 2(3uVx + 2PvUx — Su^Vx — GuvUx + 2aUxV2x — 2avu3x, (70) 

which is of order Mi = 3 in m and order M2 — 2mv. Using (34) and (35), one readily obtains 

luE — 2(3uv — 6u^v + 2aUxVx — 2au2xV, (71) 

and 

lyE = 2(3uv — 3u^v + 2aUxVx — 2au2xV- (72) 

Hence, using (33), 

J = Hu{x)E= I {luE + I,E)[Xu] — 
1 

(4/9 Amv — QX^u^v + AaXuxVx — AaXu2xv) dX 
— 2(5uv — 3u^v + 2aUxVx — 2au2xV, (73) 
which is J^^) in (59). One can set P — 1 at the end of the computations. 

6 Conservation Laws of Systems of PDEs with Trans- 
cendental Nonlinearities 

We now turn to the symbolic computation of conservation laws of certain classes of PDEs 
with transcendental nonlinearities. We only consider PDEs where the transcendental func- 
tions act on one dependent variable u (and not on x— derivatives of u). In contrast to the 
examples in the previous sections, the candidate density will no longer have constant unde- 
termined coefficients but functional coefficients which depend on the variable u. Furthermore, 
we consider only PDEs which have one type of nonlinearity. For example, sine, or cosine, or 
exponential terms are fine but not a mixture of these functions. 

6.1 The sine-Gordon Equation: Dilation Invariance and Conser- 
vation Laws 

The sine-Gordon (sG) equation appears in the literature [17, 69] in two different ways: 

• In light-cone coordinates the sG equation, Uxt = sinu, has a mixed derivative term, which 
complicates matters. We return to this type of equation in Section 7.1. 

• The sG equation in laboratory coordinates, U2t — U2x = sinu, can be recast as 

ut + v ^0, vt + U2x + smu^ 0, (74) 
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where v{x,t) is an auxiliary variable. System (74) is amenable to our approach, subject to 
modifications to accommodate the transcendental nonlinearity. 

The sG equation describes the propagation of crystal dislocations, superconductivity in 
a Josephson junction, and ultra-short optical pulse propagation in a resonant medium [69]. 
In mathematics, the sG equation is long known in the differential geometry of surfaces of 
constant negative Gaussian curvature [30, 80]. 

The sine-Gordon equation (74) is not uniform in rank unless we replace it by 

Ut + V — 0, Vt + U2x + Oisinu — 0, (75) 
where a is a real parameter with weight. Indeed, substituting the Maclaurin series, sinw = 

3 5 

u — ^ + ^ — • • ■ , and requiring uniformity in rank yields 

W{u) + W{d/dt) = W{v), 

W{v) + W{d/dt) = W{u)+2 = W{a)+W{u) = W{a) + 3W{u) = W{a) + 5W{u) = --- .(76) 

This forces us to set W{u) = and W{a) = 2. Consequently, (75) is scaling invariant under 
the symmetry 

{x, t, u, V, a) — > {X~^x, X~^t, X^u, X^v, A^a), (77) 

corresponding to W(d/dx) = W(d/dt) = 1, W(u) = 0, W(v) = 1, W(a) = 2. The first and 
second equations in (75) are uniform of ranks 1 and 2, respectively. 

The first few (of infinitely many) density-fiux pairs [8, 29] for the sG equation (75) are 

p(^) = 2acosu + v^ + ul, J^^^^2vux, (78) 

= 2vux, J^^^ ^ -2a cos u + v'^ + ul, (79) 

p^^^ — Gavux cosu + v^Ux + vu^ — 8vxU2x, (80) 

p(^) = 2q;^ cos^ u - 2o? sin^ u \.av^ cos u^v'^ ^ 20Q;ti^ cos u ^v^u\ 

Wx - - le^^L, (81) 

J^^^ and J^^^ are not shown due to length. Again, all densities and fluxes are uniform in rank 
(before a is set equal to 1). 

6.2 Computation of a Conservation Law for the sine- Gordon Sys- 
tem 

We show how to compute densities p^^^ and both of rank 2, and their associated fluxes 

j(^) and J*^^\ The candidate density will no longer have constant undetermined coefficients 
Cj but functional coefficients hi{u) which depend on the variable with weight zero [8]. To 
avoid having to solve PDEs, we tacitly assume that there is only one dependent variable with 
weight zero. As before, the algorithm proceeds in three steps: 

Step 1: Construct the form of the density 

Augment the set of dependent variables with a (with non-zero weight) and replace u by Ux 
(since W{u) — 0). Hence, V = {a^Uxiv}. Compute TZ — {a,v^ ,v'^ ,U2x-iVUx.iu'^ and remove 
divergences and equivalent terms to get S — {a, v^, u^, vUx}- The candidate density 

p — ah\{u) -\- h2{u)v^ + h3{u)u1 + hi{u)vux, (82) 
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with undetermined functional coefficients hi{u). 
Step 2: Compute the functions hi{u) 
Compute 

dp dp dp 
Dtp = -^\ut + ^ — D^ut + -^\vt 

OU OUx ov 

= {ah'i + + ul.h'^ + vuxh'^ut + {^u^hz + vhi)utx + (2't;/i2 + Uxh4)vt, (83) 
wfiere /i- denotes After replacing and Vt from (75), £^ = — D^p becomes 

£^ = {ah'i + v^h^ + tt^/^s + vuxh'^)v + {2uxhz + vh^jvx + {Ivh^ + Ma;/i4)(Q; sin m + M22;)- (84) 

£^ must be exact. Therefore, require that C^^^E = and df).E = 0. Set the coefficients of 
like terms equal to zero to get a mixed hnear system of algebraic equations and ODEs: 

h2{u) - h3{u) = 0, h'2{u)=0, h'^{u) = 0, /i1(m)=0, /i2(«)=0, 

h'iiu) = 0, 2h'2{u) - h'^{u) = 0, 2h2{u) - hl{u) = 0, (85) 

h[{u) + 2/i2(m) sin?7, = 0, h'[{u) + 2/^2 (m) sinw + 2/i2(m) cosm = 0. 

Solve the system [8] and substitute the solution 

hi{u) — 2ci cosu + C3, h2{u) — h^lu) — ci, /i4(w) = C2, (86) 

(with arbitrary constants q) into (82) to obtain 

p = ci {2a cos + t;^ + u^.) + C2t'Ma; + C3Q;. (87) 

Step 3: Compute the flux J 

Compute the flux corresponding to p in (87). Substitute (86) into (84), to get 

E — ci{2uxVx + 2vu2x) + C2{a.Ux svau + vvx + UxU2x)- (88) 

Since E—DxJ, one must integrate to obtain J. Using (26) and (35) one gets IuE — 2civUx + 
C2{ausmu + u^) and IvE — 2civUx + C2V^. Using (33), 



{IuE + I,E)[Xu] — 







{^iciXvUx + C2(Q;Msin(AM) + Xv^ + Aw^)) dX 



= Ci(2fMa,) + C2 ^-acosM + + ^M^j . (89) 

Finally, split density (87) and flux (89) into independent pieces (for ci and C2) to get 

p^^^ ^ 2acosu + v'^ + ul, J'^^^^2vux, (90) 

1 1 

= yy^^ j(2) = -a COS U + -v"^ + -M^. (91) 

For in (88), J in (89) can easily be computed by hand [8]. However, the computation of 
fluxes corresponding to densities of ranks > 2 is cumbersome and requires integration with 
the homotopy operator. 
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7 Conservation Laws of Scalar Equations with Trans- 
cendental and Mixed Derivative Terms 



Our method to compute densities and fluxes of scalar equations with transcendental terms 
and a mixed derivative term (i.e. Uxt) is an adaptation of the technique shown in Section 5. 
We only consider single PDEs with one type of transcendental nonhnearity. Since we are no 
longer dealing with evolution equations, densities and fluxes could dependent on ut,U2t, etc. 
We do not cover such cases; instead, we refer the reader to [101]. 

7.1 The sine-Gordon Equation in Light-Cone Coordinates 

In light-cone coordinates (or characteristic coordinates) the sG equation, 

Uxt = sin u, (92) 

has a mixed derivative as well as a transcendental term. A change of variables, ^ = u^,"^ — 
— 1 + cos allows one to replace (92) by 

- $ - $^ = 0, 2^ + ^2 + = 0, (93) 

without transcendental terms. Unfortunately neither (92) nor (93) can be written as a system 
of evolution equations. As shown in Section 6.1, to deal with the transcendental nonhnearity, 
which imposes W{u) = 0, one has to replace (92) by 

Uxt = Oi sin u, (94) 

where a is an auxiliary parameter with weight. Indeed, (94) is dilation invariant under the 
scaling symmetry 

(x, t, -u, a) (A^^x, \^^t, \^u, A^a), (95) 

corresponding to W{d/dx) = W{d/dt) = 1, W{u) = 0, and W{a) = 2. The density-flux pairs 
[8, 29] of ranks 2,4,6, and 8 (which are independent of Ut,U2t, etc.), are 

p(^) = ul, = 2a cos M, (96) 

p(2) = u^-Aul^ J^^^ = Aaul cos u, (97) 

p^^^ — — 20u1ul^ + 8ul^, J^^^ = 2a{3u'^ cos u + 8ulu2x sin u — 4:ul^ cos u), (98) 
^ 5ul - 280utul^ ~ lUu^^ + 22Aulul^ - 6Aul^, (99) 
j(^) = 8a (5m^ cos u + A0u'^U2x sin u + 2Qu^ul^ cos u + I6M23, ~ ISw^Wsa; cos u 

—4:8uxU2xU3x sin u + 8ul^ cos uj . (100) 

There are infinitely many density-flux pairs (all of even rank). Since Uxt = {ux)t, one can view 
(94) as an evolution equation in a new variable, U = Ux, and construct densities as linear 
combinations with constant coefficients of monomials in U and its x— derivatives. As before, 
each monomial has a (pre-selected) rank. To accommodate the transcendental term(s) one 
might be incorrectly tempted to linearly combine such monomials with functional coefficients 
hi{u) instead of constant coefficients q. For example, however, for rank R = 2, one should take 
p = ciul instead of p = hi{u)ul, because the latter would lead to Dtp = h[utul + 2hiUxU2x 
and Ut cannot be replaced from (94). 
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7.2 Examples of Equations with Transcendental Nonlinearities 



In this section we consider additional PDEs of the form Uxt — fiu), where f{u) has trans- 
cendental terms. Using the Painleve integrability test, researchers [12] have concluded that 
the only PDEs of that type that are completely integrable are equivalent to one of the standard 
forms of the nonlinear Klein-Gordon equation [2, 12]. These standard forms (in light-cone 
coordinates) include the sine-Gordon equation, u^t — sinu, discussed in Section 6.1, the 
sinh-Gordon equation u^t — sinhii, the Liouville equation u^t — e", and the double Liouville 
equations, Uxt = e" ± c~^". The latter is also referred to in the literature as the Tzetzeica and 
Mikhailov equations. For each of these equations one can compute conservation laws with 
the method discussed in Section 7.1. Alternatively, if these equations were transformed into 
laboratory coordinates, one would apply the method of Section 6.2. The multiple sine-Gordon 
equations, e.g. u^t — sin it -|- sin 2u, have only a finite number of conservation laws and are 
not completely integrable, as supported by other evidence [2]. 

The sinh-Gordon equation, u^t = sinh-u, arises as a special case of the Toda lattice 
discussed in Section 11.1. It also describes the dynamics of strings in constant curvature 
space-times [70]. In thermodynamics, the sinh-Gordon equation can be used to calculate 
partition and correlation functions, and thus support Langevin simulations [64]. In Table 1, 
we show a few density-fiux pairs for the sinh-Gordon equation in light-cone coordinates, 
Uxt = asinhw. As with the sG equation (92), W{d/dx) = W{d/dt) = 1, W{u) = 0, and 
(a) = 2. The ranks in the first column of Table 1 correspond to the ranks of the densities, 
which are polynomial in U = and its x— derivatives. The sinh-Gordon equation has 
infinitely many conservation laws and is known to be completely integrable [2]. 



Table 1: Conservation Laws of the sinh-Gordon equation, u^t — asinhu 



Rank 


Density (p) 


Flux (J) 


2 


ul 


—2acosh.u 


4 


ut + 4«2^ 


—Aau^ cosh u 


6 


ul + 20ulul^ + 8ul^ 


—2a[{3u'^ + 4:U2x) cosh-u + 8u'^U2x sinh-u] 


8 


5ul + 280ulul^ - ll2ul^ 
+224«2«2^ + 64< 


-8a [{5ul - 20ulul^ + IGulusx + 8^3^) cosh-u 
+ {A0u^U2x - 16mL + 4:8uxU2xU3x) sinhw] . 



The Liouville equation, Uxt = Q^, plays an important role in modern field theory [68], e.g. 
in the theory of strings, where the quantum Liouville field appears as a conformal anomaly 
[65]. The first few (of infinitely many) density-fiux pairs for the Liouville equation in light- 
cone coordinates, Uxt = ac", are given in Table 2. As before, W{d/dx) = W{d/dt) = 1, 
W{u) = 0, and W{a) = 2. The ranks in the table refer to the ranks of the densities. Dodd 
and BuUough [29] have shown that the Liouville equation has no densities of ranks 3, 5, and 
7. As shown in Table 2, there are two densities of rank 6, and three densities of rank 8. Our 
results agree with those in [29] , where one can also find the unique density of rank 9 and the 
four independent densities of rank 10. 
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Table 2: Conservation Laws of the Liouville equation, Uxt — oie 



R 


Density (p) 


Flux (J) 


2 




-2Q;e" 


4 




— 4Q;-u^e" 


6 


c,{ul -20ul-12ul) 


—a [Qci{ul — 4:ulu2x — S-U^a;) 

+C2U2x['iU^ + U2x)\ e 


8 


ci(-u^ - 56ulul^ - 168m^m|j. 
-672w2.wL - U4ul) 

+18m2xmL + 4mL) 
+C3{uix + 3m^mL + 15'U2,c'"L + 3mL) 


— Q! [8ci(m^ — Qu^U2x + 3m^m|j. — 20-U2X 

+C2{2ulu2x - ulul^ + 4?i^^ + 8m^-u3^ 
+24m^^X2xM3x- + 4?74,) 

+C3(4m|^ + e^^tiaj; + lSuxU2xUzx + 3M33.)] e" 



The double Liouville equations, 

^,i = e"±e-2", (101) 

arise in the field of "laser-induced vibrational predesorption of molecules physisorbed on 
insulating substrates." More precisely, (101) is used to investigate the dynamics of energy 
flow of excited admolecules on insulating substrates [82] . Double Liouville equations are also 
relevant in studies of global properties of scalar-vacuum configurations in general relativity 
and similarly systems in alternative theories of gravity [22]. 

In Table 3, we show some density-fiux pairs of u^t = a(c" — e~^"). There are no density- 
flux pairs for ranks 4 and 10. We computed a density-flux pair for rank 12 (not shown due 
to length). The results for (101) with the plus sign are similar. 

Part II: Nonlinear Differential-Difference Equations 

In the second part of this chapter we discuss two distinct methods to construct conservation 

laws of nonlinear DDEs. The first method follows closely the technique for PDEs discussed in 
Part I. It is quite effective for certain classes of DDEs, including the Kac-van Moerbeke and 
Toda lattices, but far less effective for more complicated lattices, such as the Bogoyavlenskii 
and the Gardner lattices. The latter examples are treated with a new method based on a 
leading order analysis proposed by Hickman [55]. 

8 Nonlinear DDEs and Conservation Laws 

We consider autonomous nonlinear systems of DDEs of the form 
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Table 3: Conservation Laws of the double Liouville equation, u^t = a(e" — e 



Rank 


Density (p) 


Flux (J) 


2 


ul 


-Q;(2e" + e-2«) 


4 






6 


ul + 15ulul^ - 5ul^ + 3m|^ 


-3a [{2ui + 2m2m2^ + m2Jc" 


8 


ul + A2ulul^ - Uulul^ - lu\^ 
+21«x«L - 21m2x'«L + 3mL 


-a [(8m^ + 36m^-U2x - ISm^mL - 2Qul^ 

+Quluzx + 18Ma;M2xM3x + 3Mi^)e" 

+ (4«^ - 72<M2x - l^ulul, - Aul, 

+48M^M3a; - 72Ua:U2xU3x + S^De"^"] 


10 







where u„ and F are vector-valued functions with N components. We only consider DDEs 
with one discrete variable, denoted by integer n. In many applications, n comes from a 
discretization of a space variable. The dot stands for differentiation with respect to the 
continuous variable which frequently is time, t. We assume that F is polynomial with constant 
coefficients, although this restriction can be waived for the method presented in Section 12. 
No restrictions are imposed on the degree of nonlinearity of F. If parameters arc present in 
(102), they will be denoted by lower-case Greek letters. 

F depends on Un and a finite number of forward and backward shifts of u„. We identify 
/ with the furthest negative shift of any variable in the system, and m with the furthest 
positive shift of any variable in the system. No restrictions are imposed on the integers / and 
m, which measure the degree of non-locality in (102). 

By analogy with D^;, we define the shift operator D by Du„ = u^+i. The operator D is 
often called the up-shift operator or forward- or right-shift operator. Its inverse, D^^, is the 
down-shift operator oi backward- or left-shift operator, D^^u„ = u„_i. The action of the shift 
operators is extended to functions by acting on their arguments. For example, 

DF(u„_/, u„_i, u„, u„+i, Un+m) — F(Du„_;, Du„_i, Du„, Du„+i, Du„+^) 

= F(u„_;+i, Un, U„+i, U„+2, Un+TO+l)- (103) 

Following [57], we generate (102) from 

lio = F(u_/, u_/+i, u_i, uo, Ul, u^_i, u„), (104) 

where u„ = D"^Uo = D"^F. To further simplify the notation, we denote the zero-shifted 
dependent variable, Uq, by u. Shifts of u are generated by repeated application of D. For 
instance, uj. = D^u. The identity operator is denoted by I, where D°u = lu = u. 
A conservation law of (104), 

Dtp + AJ = 0, (105) 
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links a conserved density p to an associated flux J, where both are scalar functions depending 
on u and its shifts. In (105), which holds on solutions of (104), is the total derivative with 
respect to time and A = D — I is the forward difference operator. 

For readability (in particular, in the examples), the components of u will be denoted by 
u,v,w, etc. In what follows we consider only autonomous functions, i.e. F,p, and J do not 
explicitly depend on t. 

The time derivatives arc defined in a similar way as in the continuous case, see (5) and 
(31). We show the discrete analog of (31). For a density p{up, Wp+i, • • • ,Ug, v^, Vr+i, ■ ■ ■ ,Vs), 
involving two dependent variables (m, v) and their shifts, the time derivative is computed as 




(106) 



since D and d/dt commute. Obviously, the difference operator extends to functions. For 
example, AJ = DJ — Jfora flux, J. 

A density is trivial if there exists a function ip so that p — Aip. Similar to the continuous 
case, we say that two densities, p^^^ and are equivalent if and only if p^^^ + kp^"^^ — Atp, 
for some ip and some non-zero scalar k. 

It is paramount that the density is free of equivalent terms for if such terms were present, 
they could be moved into the flux J. Instead of working with different densities, we will use 
the equivalence of monomial terms ti in the same density (of a flxed rank). Compositions of 
D or deflne an equivalence relation (=) on monomial terms. Simply stated, all shifted 
terms are equivalent, e.g. u^iVi = UV2 = U2V4 = u^s''^-i since 

U-lVi — UV2 — A. (u^iVi) — U2V4 — A {U1V3 + UV2 + U-1V1) — W_3V_i-|-A (u_2l' + Ii_3i'_i). (107) 

This equivalence relation holds for any function of the dependent variables, but for the con- 
struction of conserved densities we will apply it only to monomials. 

In the algorithm used in Sections 9.2, 11.2, and 11.4, we will use the following equivalence 
criterion: two monomial terms, ti and t2, are equivalent, ti = t2, if and only if ti = t2 
for some integer r. Obviously, if ti = t2, then ^1=^2 + A J for some polynomial J, which 
depends on u and its shifts. For example, U-2U = U-iUi because U-2U = D^^U-iUi. Hence, 

U-2U — U-iUi + [—U-iUi + U-2U\ — U-iUi + A J, with J = —U-2U. 

For efficiency we need a criterion to choose a unique representative from each equiva- 
lence class. There are a number of ways to do this. We define the canonical representative 
as that member that has (i) no negative shifts and (ii) a non-trivial dependence on the lo- 
cal (that is, zero-shifted) variable. For example, UU2 is the canonical representative of the 
class {• • • , U-2U, U-iUi, UU2, U1U3, • • • }. In the case of e.g. two variables {u and v), U2V is the 
canonical representative of the class {• • ■ , m_iW_3, mw_2, Wi't'-i, ?'2''- W3''i- • • • }• 

Alternatively, one could choose a variable ordering and then choose the member that 
depends on the zero-shifted variable of lowest lexicographical order. The code in [48] uses 
lexicographical ordering of the variables, i.e. u ^ v ^ w, etc. Thus, uv-2 (instead of U2v) is 
chosen as the canonical representative of {• • • , U-iV-s, uv-2, uiV-i, U2V, u^vi, ■ ■ ■}■ 
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It is easy to show [55] that if p is a density then D^p is also a density. Hence, using 
an appropriate "up-shift" all negative shifts in a density can been removed. Thus, without 
loss of generality, we may assume that a density that depends on q shifts has canonical form 

P(U,U1,--- ,Ug). 

9 The Method of Undetermined Coefficients for DDEs 

In this section we show how polynomial conservation laws can be computed for a scalar DDE, 

il = F{u_i, It-i+i, • • • , It, • • • , Um-l, Um)- (108) 

The Kac-van Moerbeke example is used to illustrate the steps. 

9.1 A Classic Example: The Kac-van Moerbeke Lattice 

The Kac-van Moerbeke (KvM) lattice [60, 62] , also known as the Volterra lattice, 

Un^Un{Un+l-Un-l), (109) 

arises in the study of Langmuir oscillations in plasmas, population dynamics, etc. Eq. (109) 
appears in the literature in various forms, including i?„ = l(exp(— i?n,_i) — exp(— i?„+i)), and 
Wn = Wn{w^_^i — w^.i), which relate to (109) by simple transformations [92]. We continue 
with (109) and, adhering to the simplified notation, write it as 

u = u{ui — U-i), (110) 

or, with the conventions adopted above, ii — u{Du — D''^u). 

Lattice (110) is invariant under the scaling symmetry (t.u) (A~^t, \u). Hence, u corre- 
sponds to one derivative with respect to t, i.e. -u ~ ^. In analogy to the continuous case, we 
define the weight W of a. variable as the exponent of A that multiplies the variable [39, 40]. 
We assume that shifts of a variable have the same weights, that is, W{u_i) = W{u) = W{ui). 
Weights of dependent variables are nonnegative and rational. The rank of a monomial equals 
the total weight of the monomial. An expression (or equation) is uniform in rank if all its 
monomial terms have equal rank. 

Applied to (110), W{d/dt) = W{Dt) = 1 and W{u) = 1. Conversely, the scaling symmetry 
can be computed with hnear algebra as follows. Setting W{d/dt) ~ 1 and requiring that (110) 
is uniform in rank yields W{u) — 2W{u). Thus, W{u) — 1, which agrees with the scahng 
symmetry. 

Many integrable nonlinear DDEs are scaling invariant. If not, they can be made so by 
extending the set of dependent variables with parameters with weights. Examples of such 
cases are given in Sections 11.3 and 13.2. 

The KvM lattice has infinitely many polynomial density-flux pairs. We give the conserved 
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densities of rank i? < 4 with associated fluxes (J*^^^ is omitted due to lengtfi): 

p(^) = u, J^^^ = -uu_u (111) 

p'^^^ = l-u^ + uui, J^^^ = -{u-iu^ -u-iuui), (112) 

p(^) = -V? ^ uu\{u ^ U\ ^ U2) ^ J^^^ = — (it_irt^ + 2rt_irt^iti + it_iitit^ + Xi_iwiiiii2), (113) 
3 

1 3 

P^'^ — -ii^ + H^Hi + -M^M^ + KMi(Mi + M2) + 1*1*11*2(1* + 1*1 + 1*2 + 1*3)- (114) 

In addition to infinitely many polynomial conserved densities, (110) has a non-polynomial 
density, p^^ — \nu with flux J^^ — —{u+u^]). We discuss the computation of non-polynomial 
densities in Section 12. 

9.2 The Method of Undetermined Coefficients AppHed to a Scalar 
Nonhnear DDE 

We outhne how densities and fluxes can be constructed for a scalar DDE (108). Using (110) 
as an example, we compute p^^ of rank = 3 and associated flux of rank = 4, both 
listed in (113). 

• Select the rank R of p. Start from the set V of dependent variables (including parameters 
with weight, when applicable), and form a set M. of all non-constant monomials of rank R or 
less (without shifts). For each monomial in M. introduce the right number of derivatives 
to adjust the rank of that term. Using the DDE, evaluate the derivatives, strip off the 
numerical coefficients, and gather the resulting terms in a set 7^. For the KvM lattice (110), 
V = {u\ and M. — Since and u have ranks 3,2, and 1, respectively, one 

computes 



= w^, —r- — 2uu — 2u^{ui — U-i) — 2u^ui — 2u-iu^, (US) 



and 



d^o ' dt 



d^u dii d(u(ui — u_i)) 

= u{Ui — U-i) + u{ui — M_i 



dt^ dt dt 

— u{ui — U-iY + u{ui{u2 — u) — U-i{u — U-2)) 

— uu\ — 2u^iUUi + utiU + UU1U2 — 1*^1*1 — u^iU^ + U-2U-1U, (116) 

where (110) (and its shifts) has been used to remove the time derivatives. Build TZ using the 
terms from the right hand sides of the equations in (115) and ignoring numerical coefficients, 

TZ= {u^, U^Ui, U-iU^, Uul, U-iUUi, ut.iU, UU1U2, w_2i*-ii*}. (11''') 

• Identify the elements in TZ that belong to the same equivalence classes, replace them by their 
canonical representatives, and remove all duplicates. Call the resulting set S, which has the 
building blocks of a candidate density. Continuing with (117), ■u_2i*-ii* = i*-ii*i*i = uuiU2- 
Likewise, u^iu^ = uul and u^.^u = u^Ui. Thus, S — {u^,u'^Ui,uul,uuiU2}- 
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• Form an arbitrary linear combination of the elements in S. This is the candidate p. Con- 
tinuing with the example, 

p — CiU^ + C2 U^Ui + C3 Uul + C4 UUiU2- (US) 

• Compute 




where q is the highest shift in p. Using (118) where q — 2, 

*^ \du dui du2 J 

— ((3ciit^ + 2c2UUi + csul + 041*11*2) I + (c2ii^ + Scsuui + 04x^1*2) D-|-C4ititiD^) ii. (120) 

• Evaluate Dtp on the DDE (108) by replacing u by F. Call the result E. In (110), F = 
u{ui — U-i). The evaluated form of (120) is 

E — {SciU^ + 2c2UUi + Czu\ + CiUiU2)u{ui — W_i) + {C2U^ + 2c^UUi + C4UU2)ui{u2 — u) 
-\-{c4UUi)u2{u^ - Ui) 
= (3ci — C2)U U\ — 'iC\U-\U + 2(C2 — C-i)U Ui — 2c2U-iU Ui + C3UU1 — C3U-1UU1 

— C4U-1UU1U2 + (C2 — Ci)U^UiU2 + 2C2,UU\U2 + C^UUiU^ + C4WWiIi2li3- (121) 

• Set J = 0. Transform £^ into its canonical form. In doing so modify J so that E + AJ 
remains unchanged. For example in (121), replace — 3ciU_iU'^ in E by —'iciuu\ and add 
—'iciUiV? to J since liwf — [uu\ — u^iu^]. Do the same for all the other terms which are not 
in canonical form. After grouping like terms, (121) becomes 

E = (3ci - C2)u^ui + (c3 - 'ici)uu\ + 2(C2 - C3)u^u\ + 2(C3 - C2)uu\u2 

+ (C4 - C3)UUIU2 + (C2 - C4)m^MiM2, (122) 

with 

J = — (3cili_ili^ + 2c2li-lli^Wl + CzU-iUu\ + C4li_ililiili2)- (123) 

• £■ is now the obstruction to p being a density. Set E — and solve for the undetermined 
coefficients Cj. Thus, 

3ci - C2 = 0, 3ci - ca = 0, C2 - C3 = 0, C3 - C4 = 0, C2 - C4 = 0, (124) 

which yields 02 = 03 = 04 = 3ci, where Ci is arbitrary. 

• Substitute the solution for the q into the candidates for p and J to obtain the final density 
and associated flux (up to a common arbitrary factor which can be set to 1 or any other 
nonzero value). For the example, setting Ci = | and substituting C2 = C3 = C4 = 1 into (118) 
and (123) yields p(^) and J^^) as given in (113). 
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10 Discrete Euler and Homotopy Operators 



For simplicity, we will consider the case of only one discrete (dependent) variable u. First, we 
remove negative shifts from E. Thus, E — D'^E where I corresponds to the lowest shift in E. 

The discrete variational derivative (discrete Euler operator) 

An expression E is exact if and only if it is a total difference. The following exactness test 
is well-known [10, 57]: A necessary and sufficient condition for a function E, with positive 
shifts, to be exact is that Cu^ E = 0, where jCu^ is the discrete variational derivative (discrete 
Euler operator of order zero) [10] defined by 

C(^)E = |- (f^D"M£; = |-(l + D-^ + D-2 + D-3 + ... + 

^ \k=o J ^ 

,dE ^ ,dE ^ ,dE ^ .dE m dE , , 

= 1^ + D-'— + D-2— + D-3— + . . . + D-'"-— , 125 

OU OUi OU2 OU3 aum 

where m is highest shift occurring in E. 

Application. We return to (121) where 1 = 1. Therefore, 

E = DE = (3Ci — C2)u^U2 — 3CiUU\ + 2(C2 — C3)U^U2 — 2C2UU]U2 + ~ C3UU1U2 

— C4UU1U2U3 + (C2 — C4)ulu2Us + 2c3Uiulu3 + C4U1U2UI + C4U1U2U3U4, (126) 

which is free of negative shifts. Applying (125) to (126), where m = 4, gives 

40)^ = 1- (I + D-^ + + + D-^) ^ 

= 3(3ci — C2)u'^Ui + 3(c3 — 3ci)-u_i-u^ + 2(c2 — Ci)uuiU2 + 4(c2 — c^)uu\ 

+4(C3 - C2)u_iUUi + 2(C3 - C2)u\u2 + (C3 - 'iCi)u\ + (C4 - C^)u_iu\ 

+ (C4 - C-i)uiul + (3Ci - C2)u\ + (C2 - C4)u\ui + 4(c2 - C3)m?_iM 

+2(C3 - C2)U_2U\ + 2(C4 - C3)m_2M-iM + (C2 - C4)u\u_i, (127) 

which, as expected, vanishes identically when Ci = |, C2 = C3 = C4 = 1. 

Due to the large amount of terms generated by the Euler operator, this method for finding 
the undetermined coefficients is much less efficient than the "splitting and shifting" approach 
illustrated on the same example in Section 9.2. 

The discrete homotopy operator 

As in the continuous case, the discrete homotopy operator reduces the inversion of the dif- 
ference operator, A = D — I, to a problem of single- variable calculus. Indeed, assuming that 
E is exact and free of negative shifts, the flux J = A~^£' can be computed without "summa- 
tion by parts." Instead, one computes a single integral with respect to an auxiliary variable 
denoted by A (not to be confused with A in Section 9.1). 

Consider an exact expression E (of one variable m), free of negative shifts, and with highest 
shift m. The discrete homotopy operator is deflned [61, 71, 72] by 

nuE= / {IuE)[Xu]—, (128) 
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with 

m 

4^ = E I E ) D-^^ = E I E D-M (129) 




A:=l 

where {IuE)[Xu] means that in I^E one replaces u Xu, Ui — > A^i, U2 \u2, etc. The 
formulas in (129) are equivalent to the one in [52], which in turn is equivalent to the formula 
in terms of discrete higher Euler operators [51, 54]. Given an exact function E without 
negative shifts one has J = A''^E — HuE. A proof can be found in [61, 72]. 

Application. Upon substitution of Ci = |, C2 = C3 = C4 = 1 into (126), we obtain 

E=DE= —uu\ — 2uu\u2 + Uiu\ — uuiu\ — UU1U2U3 + 2uiu\u^ + UiU2u\ + (130) 
where the highest shift is m = 4. Using (129), 

-1 , n-2 , r.-3\„. , /n-l , nv-2 , n>-3 , n-4x„. ^-^ 



' du2, du4 

I/O QO 00 \ 

= D~ ui{—3uui — AuuiU2 + U2 — UU2 — UU2US + 2U2US + U2U^ + 1*21*31*4) 

+ (D~^ + D~'^)u2{ — 2uu\ + 2>Uiu\ — 2uUiU2 — UUiUs + 4l*iU2l*3 + Ulul + U1U3U4) 

+ (D~"^ + D^^ + D^'^)us{—uuiU2 + 2uiu\ + 2u]U2Uz + 1*11*21*4) 
+ (D-^ + + + D-^)i*4(i*iM2i*3) 

= A{uu\ -\- 2uu\u2 -\- UUiu\ -\- UUiU2U'i) ^ (131) 

which has the correct terms of J^^^ but incorrect coefficients. Finally, using (128) 

~ fix 

J = nu{-E) = - j^{IuE)[Xu] — 

— —4 / (^uu\ + 2uu\u2 + uuiu\ + uuiU2U'i) X^ dX 

Jo 

— —{uul + 2uulu2 + uuiul + uuiU2U2)- (132) 
Recall that E — DE. Hence, 

J = D~^J — —{u-iu^ + 2u-iu'^ui + U-iuu1 + U-1UU1U2), (133) 

which corresponds to J^^^ in (11). 

The homotopy method is computationally inefficient. Even for a simple example, like 
(131), the integrand has a large number of terms, most of which eventually cancel. To 
compute the flux we recommend the "splitting and shifting" approach which was illustrated 
(on the same example) in Section 9.2. 

The generalization of the exactness test to an expression E with multiple dependent 
variables {u,v,- • •) is straightforward. For example, an expression E of discrete variables u, v 
and their forward shifts will be exact if and only if Cu^E — (Cu^ E , cij'^ E) = (0,0), where 
C^^ is defined analogously to (125). Similar to the continuous case, the homotopy operator 
formulas (128) and (129) straightforwardly generalize to multiple dependent variables. The 
reader is referred to [51, 52, 54] for details. 
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11 Conservation Laws of Nonlinear Systems of DDEs 

We use the method discussed in Section 9.2 to compute conservation laws for the Toda and 
Ablowitz-Ladik lattices. Using the latter lattice, we illustrate a "divide and conquer" strategy, 
based on multiple scales, which allows on to circumvent difficulties in the computation of 
densities of DDEs that are not dilation invariant. 

11.1 The Toda Lattice 

One of the earliest and most famous examples of completely integrable DDEs is the Toda 
lattice [93, 94], 

ijn = CXp {yn-1 - Vn) " exp - Vn+l), (134) 

where y„ is the displacement from equilibrium of the nth particle with unit mass under 
an exponential decaying interaction force between nearest neighbors. With the change of 
variables, w„ = f„ = exp (y„ — yn+i), lattice (134) can be written in algebraic form 

iln^Vn-l-Vn, Vn ^ Vn{Un - Un+l) ■ (135) 

Adhering to the simplified notation, we continue with 

ii — V-i — v, v — v{u — ui). (136) 

As before, we set W{d/dt) — 1, assign unknown weights, W{u),W{v), to the dependent 
variables, and require that each equation in (136) is uniform in rank. This yields 

W{u) + l = W{v), W{v) + l = W{u) + W{v). (137) 

The solution W{u) — 1, W{v) — 2 reveals that (136) is invariant under the scaling symmetry 

{t, u, v) {X-H, Xu, X'^v), (138) 

where A is an arbitrary parameter. 

The Toda lattice has infinitely many conservation laws [44]. The first two density-fiux 
pairs are easy to compute by hand. Here we give the densities of rank i? < 4 with associated 



fluxes, J(^^ being omitted due to length: 

p« = u, J^^^=v^i, (139) 

= -u'' + v, J(^)=ii^;_i, (140) 

= J^u{V-i+v), J'^^^ =U-iUV-i + v\, (141) 

= }-u^ j^u'^(^y_-^j^y'j j^uuiv + ^v'^ + VVl. (142) 
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11.2 Computation of a Conservation Law of the Toda Lattice 

As an example, we compute density p^^^ (of rank R — ?,) and associated flux J^^^ (of rank 4) 
in (141). 

Step 1: Construct the form of the density 

Start from V = {u, v}, i.e. the set of dependent variables with weight. List all monomials in u 
and V of rank i? = 3 or less: M. = {m^, -u^, uv, u, v}. Next, for each monomial in Ai, introduce 
the correct number of t-derivatives so that each term has rank 3. Using (136), compute 

d°«3 3 d%v 
^ = ^' ^ = ^"' 

d-u^ . ^ dv . 

—— = 2uu — 2uv-i — zuv, — — V — uv — uiv, 
dt dt 

d'^u du d(v_i — v) 



dt^ dt dt 



u^iV-i — uv^i — uv + Uiv. (143) 



Gather the terms in the right hand sides in (143) to get 71 = {u'^jUV-ijUVjU-iV-ijUiv}. 

Identify members belonging to the same equivalence classes and replace them by their 
canonical representatives. For example, uv^i = Uiv. Adhering to lexicographical ordering, 
we will use uv^i instead of Uiv. Doing so, replace 7^ by <S = {u^,uv-i,uv}, which has the 
building blocks of the density. Linearly combine the monomials in S with undetermined 
coefficients q to get the candidate density of rank 3 : 

p — ciu^ + C2UV-1 + C3UV. (144) 

Step 2: Compute the undetermined coefficients q 

Compute Dfp and use (136) to eliminate ii and v and their shifts. Thus, 

E = (3ci - C2)u^V-i + (C3 - 3ci)u^V + (Ca - C2)V-iV + C2U-1UV-1 + C2v\ 

—CsUUiV — c^v^. (145) 

Step 3: Find the associated flux J 

Transform (145) into canonical form to obtain 

E — (3ci — C2)u1v + (ca — 3ci)u^v + (ca — C2)vvi + C2UU1V + C2V^ — csuuiv — c^v^ (146) 
with 

J = (3ci - C2)u^v-i + (ca - C2)v^iv + C2U-1UV-1 + C2v\. (147) 
Set £■ = to get the hnear system 

3ci-C2 = 0, ca-3ci = 0, C2 - Ca = 0. (148) 

Set Ci = I and substitute the solution ci = |, C2 = Ca = 1, into (144) and (147) to obtain p*^^) 
and J(3) in (141). 
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11.3 The Ablowitz-Ladik Lattice 

In [4, 5, 6], Ablowitz and Ladik derived and studied the following completely integrable 
discretization of the nonlinear Schrodinger equation: 

i iln = Un+1 - 2Un + Un-l ± ulUn{Un+l + Wn-l), (149) 

where -u* is the complex conjugate of We continue with (149) with the plus sign; the 
case with the negative sign is analogous. Instead of splitting m„ into its real and imaginary 
parts, we treat Un and Vn — u"^ as independent variables and augment (149) with its complex 
conjugate equation, 

Un = {Un+l-^Un + Un-l) +UnVn{Un+l + ^ln-l), 

Vn = -{Vn+1 - 2Vn + Vn-l) - Mn^^„(t^n+1 + Vn-l), (150) 

where i has been absorbed into a scale on t. Since Vn = "U* we have W{vn) = W{un). Neither 
equation in (150) is dilation invariant. To circumvent this problem we introduce an auxiliary 
parameter a with weight, and replace (150) by 

it = a{ui — 2u + U-i) + uv{ui + U-i), 

V = —a{vi — 2v + v_i)—uv{vi + v_i), (151) 
presented in the simplified notation. Both equations in (151) are uniform in rank provided 

W{u) + 1 = W(a) + W{u) = 2W{u) + W{v), 

W{v) + 1 = W(a) + W{v) = 2W{v) + W{u), (152) 

which holds when W{u) + W{v) = W{a) = 1. Since W{u) = W{v) we have W{u) = W{v) — 
1, and W{a) = 1, which expresses that (151) is invariant under the scaling symmetry 

{t, u, V, a) — > {X^^t, X'^u, X^v, Xa). (153) 

We give the conserved densities of (151) of ranks 2 through 4. Only the flux of rank 3 
associated to p^^^ is shown. The others are omitted due to length. 

p(^) — a{ciuv^i + C2UV1), (154) 

J^^^ — —a{ci{auv-2 — (^u^iV-i + U-1UV-2V-1) + C2{auv — au-iVi — U-iuvvi)), (155) 

p(^) — a [ci{^u^v\ + uuiv^iv + auv_2) + C2{^u^vl + UU1V1V2 + auv2)) , (156) 

p(^) — a (ci [l-u^f^i + -UMif _if (-uf _i + Uiv + U2V1) + auv_i{uv_2 + Wif_i) 

+aUv{uiV-2 + U2V-1) + a'^UV-s] + C2 [^U^vf + UUiViV2{uVi + U1V2 + U2V'i) 

+auv2{uvi + U1V2) + auv2,{uiVi + U2V2) + a^uv^ ) , (157) 

where Ci and C2 are arbitrary constants. Our results confirm those in [6]. Since our method 
is restricted to polynomial densities we cannot compute the density with a logarithmic term, 

Pn = «(tin-l + Un+l) " 2 ln(l + Xi„0) , (158) 

which corresponds [3, 6] to the Hamiltonian of (149), viz. H — —i YlnP^i- 
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11.4 Computation of a Conservation Law of the Ablowitz-Ladik 
Lattice 

To make (150) scaling invariant we had to introduce an auxiliary parameter a. This compli- 
cates matters in two ways as we will show in this section. First, to compute rather simple 
conserved densities, like p = uv_i and p = uvi, we will have to select rank i? = 2 for which the 
candidate density has twenty terms. However, eighteen of these terms eventually drop out. 
Second, conserved densities corresponding to lower ranks might appear in the result. These 
lower-rank densities are easy to recognize for they are multiplied with arbitrary coefficients 
Cj. Consequently, when parameters with weight are introduced, the densities corresponding 
to distinct ranks are no longer linearly independent. 

We compute p^^^ of rank R = 2 in (154) with associated flux J^^^ in (155). Note that p^^^ 
and J^^^ cannot be computed with the steps below when R — 1. 

Step 1: Construct the form of the density 

Augment the set of dependent variables with the parameter a (with non-zero weight). Hence, 
V = {u, V, a}. Construct 

Ai = {u, V, u , uv^ V , ati, u , av, u v,uv ,v ,au ,u , auv, u v,av ,u v ,uv ,v j, (159) 

which contains all non-constant monomials of (chosen) rank 2 or less (without shifts). As 
with the previous examples, for each element in Ai add the right number of t— derivatives. 
Use (151) to evaluate the t— derivatives, gather the terms in the right hand sides, introduce 
the canonical representatives (based on lexicographical ordering), and remove duplicates to 
get 

S = {au'^,u^,auui,auv_i,auv,u^v,u'^uiv,u'^v_iv,av'^,u^v'^, 

uuiv ,uv-iv ,uv ,v ,auvi,uUiVi,avvi,u vvi,uv vi,uuiVi}. (160) 

Linearly combine the monomials in S with undetermined coefficients Cj to get the candidate 
density of rank 2: 

p — Ci mi? -|- C2 -|- C3 auu\ -|- C4 cmv^x -|- C5 auv -|- -|- C7 u^uiv 

+Cs U^V-iV -\- Cg aV^ -\- Clo U^V^ -\- Cii UUiV^ + C12 UV-iV^ -\- Ci3 UV^ -\- Ci4 

-I-C15 aUVi -\- C16 UU]Vi -\- C17 aVVi -\- Ci8 U VVi + +Cig UV Vi + C20 UUiV^ . (161) 

Step 2: Compute the undetermined coefRcients q 

The computations proceed as in the examples in Sections 9.2 and 11.2. Thus, compute 
DjP and use (151) to eliminate ii and v and their shifts. Next, bring the expression E into 
canonical form to obtain the linear system for the undetermined coefficients q. 

Without showing the lengthy computations, one finds that all constants q = 0, except C4 
and C15 which are arbitrary. Substitute the coefficients into (161) to get p^^^ in (154). 

Step 3: Find the associated flux J 

The associated fiux comes for free when E is transformed into canonical form. Alternatively, 
one could apply the homotopy approach for multiple dependent variables [51, 54] to compute 
the flux. In either case, one gets J^^^ in (155). 
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11.5 A "Divide and Conquer" Strategy 

It should be clear from the example in Section 11.4 that our method is not efficient if the 

densities arc not of the form (111)-(114) for the KvM lattice and (139)-(142) for the Toda 
lattice. Indeed, the densities for the AL lattice in (154)-(157) are quite different in structure. 
Therefore, in [34], Eklund presented alternate strategies to deal more efficiently with DDEs, 
in particular with those involving weighted parameters such as (151). 

A first alternative is to work with multiple scales by setting either VF(Dt) = or W{Dt) — 
1, the latter choice is what wc have used thus far. A second possibility is to leave VF(D() 
unspecified and, if needed, introduce extra parameters with weight into the DDE. 

Let us explore these ideas for the AL lattice (151), which we therefore replace by 

u — aP{ui — 2u + u^i) + (3uv{ui + u^i), 

V — —a(3{vi — 2v + V-i)—Puv{vi + V-i), (162) 

where /3 is a second auxiliary parameter with weight. Requiring uniformity in rank leads to 

W(u) + W(Dt) = W(a) + W(P) + W(u) ^W((3) + 2W(u) + W{v), 

W(v) + W(Dt) = W{a) + W((3) + W(v) ^W((3) + 2W(v) + W(u), (163) 

which are satisfied when W{a) = W{u)+W{v) and W{Dt) = W{u)+W{v)+W{l3). Therefore, 
we can set W{Dt) — a, W{u) — 6, and W{P) — c with a, 6, c rational numbers so that 
W{v) — a — b — c and W{a) = a — c are strictly positive. Thus (162) is dilation invariant 
under a three-parameter family of scaling symmetries, 

(t, u, V, a, (3) (X-% X\ X^-^-^v, X^-^a, A"/3). (164) 

Scaling symmetry (153) corresponds to the case where a = 1, 6 = |, and c = 0, more precisely, 
j3 = 1. The fact that (162) is invariant under multiple scales is advantageous. Indeed, one 
can use the invariance under one scale to construct a candidate density and, subsequently, 
use additional scale(s) to split p into smaller densities. This "divide and conquer" strategy 
drastically reduces the complexity of the computations. The use of multiple scales has proven 
to be successful in the computation of conservation laws for PDEs with more than one space 
variable [51]. 



i 


Rank 


Candidate pi 


Final pi 


Final Jj 


1 


a + 2b-c 


Ciav? + c^auui + CqW'v + cru^uiv + CiqUuIvi 








2 


Ah 


4 

C2M 








3 


2{a - c) 


c^auv-i + c^auv + csu^V-iv + c\qV?v^ 

+CiiUUiV^ + Ci^aUVi + CisU^VVi + C2qUUiv\ 


C4auv-i 
+Ci5auvi 


J^i) in (155) 


4 


3a — 2b — 3c 


cgav"^ + ci2UV-iv^ + cisuv^ + ciravvi + ciguv^vi 








5 


4{a-b- c) 


4 









Table 4: Candidate densities for the Ablowitz-Ladik lattice (162). 
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Candidate density (161) is uniform of rank 2 under (153) but can be split into smaller 
pieces, pj, using (164), even without specifying values for a,b, and c. Indeed, as shown in 
Table 4, p in (161) can be split into pi through p^ of distinct ranks under (164). Steps 2 
and 3 of the algorithm are then carried out for each of these pi, {i = 1, ■ ■ ■ , 5) separately. As 
shown in the table, all but one density lead to a trivial result. The longest density, pa, with 
8 terms, leads to p^^^ and J^^^ in (154) and (155), respectively. 

12 A New Method to Compute Conservation Laws of 
Nonlinear DDEs 

In the continuous case, the total derivative has weight one. Consequently, any density is 
bounded with respect to the order in x. For example, as shown in Section 3.2 for the KdV 
case, the candidate density p of rank 6 in (18) is of first order (after removing the second and 
fourth order terms.) Obviously, a density of rank 6 could never have leading order terms of 
fifth or higher order in x because the rank of such terms would exceed 6. 

In the discrete case, however, the shift operator D has no weight. Thus, any shift D'^u = Uk 
of a dependent variable u has the same weight as that dependent variable. Simply put, 
W{uk) — W{u) for any integer k. Consequently, using the total derivative as a tool 
to construct a (candidate) density has a major shortcoming: the density may lack terms 
involving sufficiently high shifts of the variables. As shown in Section 9.2 for the KvM 
lattice, the candidate density p of rank 3 in (118) has leading order term uuiU2, i.e. the term 
with the highest shift (2 in this example). It is a priori not excluded that p might have 
terms involving higher shifts. For example, uuiu^ has rank 3 and so do infinitely many other 
cubic terms. Note that for this example we constructed p starting from powers of u, viz. 
u^,u'^,u; and, by repeated difi^erentiation, worked our way "down" towards the leading order 
term uuiU2- In this section we outline the key features of a new method which goes in the 
opposite direction: (i) first compute the leading order term and subsequently (ii) compute 
the terms involving lower shifts. In step (ii) only the necessary terms are computed, nothing 
more, nothing less. This method is fast and powerful for it circumvents the use of the dilation 
invariance and the method of undetermined coefficients. More importantly, the new method 
is not restricted to densities and fluxes of polynomial form. 



12.1 Leading Order Analysis 

Consider a density, p, that depends on q shifts. Since D'^ p is also a density, we may, without 
loss of generality, assume that p has canonical form p(u, ui, . . . , Uq) with ^^^^ 7^ 0. In 
[55], Hickman derived necessary conditions on this leading term (which, in the system case, 
is a matrix) . First all terms in the candidate density p that contribute directly to the flux are 
removed. Rather than applying the Euler operator on the remaining terms in p, the necessary 
condition [57], 

^ = 0. (165) 

for to be a total difference is applied to obtain a system of equations for the terms that 
depend on the maximal shift, u^, in p. This system is rewritten as a matrix equation. Solutions 
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to this system will give us the form of the leading term in p. 
We apply a splitting of the identity operator, 



l = (D-l + l)D-^ = AD-^ + D-\ (166) 

to the part, p*, of the candidate p that is independent of the variables with the lowest order 
shift. The first term AD^^p* contributes to the flux while the second term D~^p* has a 
strictly lower shift than p* . Applying this split repeatedly we get 

I - (D^ - I + I) D"^ = A (D'^-^ + D'^-^ + . . . + D + I) D"^ + D-^ (167) 

where, again, the first term contributes to the flux and the remainder has strictly lower shift. 

This decomposition is repeatedly applied to terms that do not involve the lowest order 
shifted variables. Any terms that remain will involve the lowest order shifted variable. These 
terms yield the constraints on the undetermined coefficients or unknown functions in the 
density p. As shown in [55] , the result of this split is that p is a density if and only if 

-D'(F^)tDV+ED'(F|;)p (168) 

\ / j=0 3=1+1 ^ 

is a total difference, where the operator 

F^ = yF"-^ (169) 
dn ^ du'^ ^ ^ 

a 

involves a summation over the components in the system of DDEs. As before, in the examples 
we will use u,v,w, etc. to denote the dependent variables u'^. 

Now, (7 depends on the shifted variables u, . . . , Ug+L where L — max {I, m). For q > L, 
(165) gives 

D^f^^Vo'fl^V^^O, (170, 



dudugj^L ydudugdu^L/ \duL J dudug 

where T stands for transpose. The system case is treated in detail in [55]. For brevity, we 
continue with the scalar case. Let 

OF 

dU-L OUl 

then the leading term will satisfy 

with S = A + We immediately see that if / 7^ m then either A or p is zero and 

(172) has no non-trivial solutions. Let q = pL + r with p and r integers, < r < L, and 

Yl D'^^A D^'^C- (173) 

Kk=l J 



38 



Then 

p-i 



Sc = Yl D'^-^A (A D^C + C D» . (174) 

k=l 

Thus c ^ lies in the kernel of S if and only if 

A D^C + C D> = (175) 
has a non-zero solution (. Suppose (175) has two non-zero solutions, say, and (2- Then, 

I ^ 0^(1). 

So, since L ^ 0, (2 = ^Ci fo^^ some constant a. Therefore, the kernel of »S is, at most, one 
dim,ensional. This implies that a scalar DDE can have, at most, one conserved density that 
depends on q shifts for q > L. The leading term, p, will satisfy 

c, (177) 



du dUq 



which, upon integration, yields 



cduduq. (178) 



The density (if it exists) may now be computed by a "split and shift" strategy on this leading 
term. Starting with p = p. the objective is to successively compute the terms (of lower shift!) 
that must be added to p until p = 0. First, p is computed and evaluated on the DDE. 
Next, all terms are shifted so that the resulting expression depends on u (and not on lower 
shifts of u). Then the leading terms, ^, in p are isolated. Last, we solve 

Dt p^^^ = C + terms of lower shift. (179) 

If (179) has no solution then a density with q shifts does not exist. On the other hand, if (179) 
has a solution, the "correction" term p^^^ is then subtracted from p and we recompute p. 
By construction, the highest shift that occurs in the result will now be lower than before and 
we repeat the entire procedure to obtain a new correction term p*^^\ After a finite number of 
steps, we will either find an that the correction term does not exist (and so the density does 
not exist) or we will obtain p = and p will be a density. 

This algorithm has been coded [56] in Maple. Further details and worked examples will 
be presented in [58]. For now, we illustrate the algorithm with a simple example. 

12.2 A "Modified" Volterra Lattice 

Consider the DDE, 

u^u'^{u2-U-2), (180) 
which is related to the well-known modified Volterra Lattice [11]. Here L — 2 and, using 
(171), 

A=|^=-^.^ H-^ = u'. (181) 
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Thus, the condition (175) for a non-trivial density becomes C,u1 = D^^ for r = 0, 1. For 
r = 0, we have C, = and so we may choose C = 1- For r = 1, we have C,u\ = D^(^, which 
has no non-zero solutions. Thus, densities only exist for r = 0. Since q = pL = 2p, with p 
integer, we conclude that q must be even. In these cases the leading term will satisfy 



^2 /p-i 



du duq 



Y[u2k] ^UU2 ■■■ (182) 



Therefore, after scaling to remove constants, the leading term is p — u^U2 • • • Uq^2 Uq- 

For example, let us compute the density that depends on g = 4 shifts. We set p = p = 
u'^U2U4. Using (180), 

DtP — U U2U4{U2 — U-2) + 2UU2U4{U4 — U) + UU2U4^{Uq — U2) 

= uu^ul — U2U4. (183) 
The highest shift is U4 and so the leading terms are 

^ = uulul- U2^ U4. (184) 

Note that the terms in must cancel by the construction of p. Next, we note that U4 as a 
leading term must arise as either U2 = ul {U4 — u) or 

U2U — U2U {u2 — U-2) = u U2 — UU2U4. (185) 

The quadratic term must arise from (185) (as we have already determined all terms that 
involve U4, so, we cannot have a quadratic term given by U4U2). Now, we solve (179) to get 

p(i)^_i^,2^2^... . (186) 

Indeed, 

D/ 1 2 2^ 2 • 2 • 

t(— 2^ U2) = —UU2U — U U2U2 

= uu\u\ — u\ U4 + lower shift terms 

— ^ + lower shift terms, (187) 

with ^ in (184). Thus, we update p — u'^U2U4 by subtracting p'^^\ This yields, 

p — U2U4-\- ^v?' li^. (188) 

We readily verify that D^p = and, thus, p in (188) is the unique density of (180) that 
depends on 4 shifts. 

13 The Gardner Lattice 

In this section we consider the DDE described in [91], 

(1 + ah?u + (5h?u^) |^ (i'^-2 - i^-i + tii - \u'^ + ^ [u-iU-2 + u\ + u{u-i - ui) 



u 



-uj - U1U2] + ^ [Mii(ti-2 + u)- ul{u + U2)] I , 



(189) 
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which is a completely integrable discretization of the Gardner equation, 

Ut + 6auUx + SPu^Ux + — 0. (190) 

Therefore, we call (189) the Gardner lattice. Note that (190) is a combination of the KdV 
equation (/3 = 0) and the mKdV equation {a — 0). Consequently, (189) includes the com- 
pletely integrable discretizations of the KdV and mKdV equations as special cases. 



13.1 Derivation of the Gardner Lattice 

Based on work by Taha [91], we outline the derivation of the Gardner lattice from a discretized 
version [4] of the eigenvalue problem of Zakharov and Shabat. Consider the discrete system 
[1, 7], 

Vi,n+l = zVi.n + Qn{t)V2,n: V^2,n+1 = ^V^2,n + K(t)^l,n, (191) 

Vi,n = An{z)V,,n + Bn{z)V2,n, l>2,n = Cn{z)V,,n + Dn{z)V2,n, (192) 

where, in general, the coefficients An through D„ depend on the potentials, Rn and Qn- The 
eigenvalue z is assumed to be time-independent. 

Step 1: Expressing the compatibility conditions, V^^n+i — DVi,„, for i — 1,2, and equating 
the coefficients of Vi^n and V2,n, we get [7] 

1 

CnQn — Bn+lRn = z{An+l — An), BnRn — Cn+lQn = -{^n+l ~ Dn), (193) 

Z 

1 ■ 1 

Qn + DnQn ~ ^n+lQn — ~Bn+\ ~ ^Bn, Rn + AnRn — Dn+lRn = ~Cn+\ ~ zCn- (194) 

z z 
Step 2: Substituting the expansions [91], 

An=Y. '"'^n\ Dn=j^ z'W^^^\ (195) 

j=-2 j=-2 

Bn-Y. ^''-'Bi^'-'\ Cn^Yl -''-'Cj^'-'\ (196) 
i=-i i=-i 

into (193) and (194) and setting the coefficients of like powers of z equal to zero, we ob- 
tain a system of twenty equations for the eighteen unknowns functions A'^'^Di^''^ with 
j = —2,-1,0,1,2; and Bn"" ^\ Cn''^^^ with j = —1,0,1,2. The simplest equations arise 
from the coefficients of the terms in z"^^ and : 

- = 0, Ai;? - Di'^ = 0, (197) 

Qn{D^:^ - Al%) + fif = 0, i?„(A(-^) - dj) + Ct'^ = 0, (198) 
Qn{D(-'^ - 4;?) - i^i;? = 0, Rn{A(:^ - D^l) - = 0, (199) 

Step 3: We outhne the solution process which follows the strategy in [7, Section 2.2a]. 
Prom (197), we conclude that An^ and d\t^^^ are independent of n- Hence, and 
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are constants. The tilde will remind us that we are dealing with constants. 
Solving (198) and (199), we get 



(200) 
(201) 



Substituting these results into two of the equations coming from z^^, we find that An - 
A^'^^ and Dlf' — D^^^ are constants. Prom equations corresponding to z^^, z^^, we obtain 



AD(-2) = Z^i;^,) - = ^{RnQn-1 - Rn+lQn), 

AAi"^) = Ai-_,'l - 4-2) = P{QnRn+l - Qn-lRn), 

AL>(^) = D^^l - L>(^) = a(it:„g„+i - Rn-iQn), 



(202) 
(203) 
(204) 
(205) 



where a = i^^^ - D'^^^ and /3 = D^'^^ - A^-^\ Equations (202)-(205) are satisfied for 

42) ^ ^(2) _ = ^(-2) _ (206) 

4-2) ^ ^(-2) ^ ^g^_^i?^^ £,(2) ^ £,(2) + «i?„_,Q„, (207) 

where A^'^\ D^'"^^ and Jl*^"^^ £>(^) are constants. Next, we solve equations (from the terms in 
z±2) for B^n^^ and C^^^\ yielding 



^i'^ = SQn + aQn+i - aQniQ n-\-lRn ~t~ QnRn—l)-} 

— + PRn+1 — i^Rn{Rn+lQn + RnQn-l), 

— 7Qn-l + f3Qn-2 — (3Qn-l{Rn-lQn-2 + RnQn-l), 

n-1 + aRn-2 — OiRn-l{Qn-lRn-2 + QnRn-l), 



(208) 

(209) 
(210) 
(211) 



where 5 = A*^^) — Z^*^^) and 'j = '^'^ —A^ Next, we solve equations coming from 2;^^, which 
involve A74n'* and ADn \ These equations, which are similar in structure to (202)- (205), yield 

~ ^ {~Qn+lRn-l ~ QnRn-2 + QnQn+lRn~lRn + Qn-lQnRn-2Rn-l + Qn-^n-l) 

- 5Rn-iQn, (212) 

— ^ {~Rn+lQn-l — RnQn-2 + RnRn+lQn-lQn + Rn-lRnQ n-2Q n-l + -^nQn-l) 

+L'(°)-7g„_ii?„, (213) 

where A^^^ and D'-^-' are constants. All the difference equations for the unknown coefficients 
are now satisfied. We are left with the two DDEs (coming from the terms in z^), 

Qn - liQn + (1 - QnRn) -OiQn+2 " ^Qn+l + iQn-l + ^Qn-2 + OcQ n+l{Q n+2Rn+l 

-\-Qn+lRn) ~ f^Qn-l{Qn-2Rn-l + Qn-lRn) + OiQnQn+lRn-1 — PQuQu-lRn+l = 0, (214) 

Rn + kRu + (1 ~ RnQn) —PRn+2 — jRn+1 + ^Rn-1 + OiRn-2 + PRn+l{Rn+2Qn+l 

+Rn+lQn) — OlRn-l{Rn-2Qn-l + Rn-lQn) + PRnRn+lQn-1 — ^RnRfi-lQn+l — 0, (215) 
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where K = - L>(°). 

Step 4: The terms —RQn in (214) and kR^ in (215) could be removed with a suitable 
transformation. Accomplishing the same, we set k = 0. Next, we substitute 

Qn = Un, Rn = -h^ia + /3Un), (216) 

into (214) and (215). Next, we set = a and 7 = 5 to remove all constant terms. Doing so, 
(214) and (215) collapse into a single DDE, 

Un + {1 + ah^Un + (3h^u'^n)^a{Un-2 - Un+2) + S{Un-l - Un+l) + aOih^ {Un-2Un-l + 
+Un{Un-l-Un+l) -^n+l -Un+lUn+2) + pah'^ {ul_^{Un-2 + U„) -ul^^{Un + Un+2))} = 0, (217) 

Step 5: To fix the scale on t, as well as the constants a, f3, and S, we consider the limit of 
(217) for h —> 0. Using, limh^QUn{t) — u{x,t), limh^oUn{t) — Ut{x,t), and substituting 

^mun+m{t) = u{x,t) + mhux{x,t) + ^{mh)'^U2x{x,t) + ^{mh)^U3x{x,t) + . . . , (218) 

into (217), we obtain 

ut - 2h{S + 2a)ua: - 2h^{S + 8a) {auu^ + /3u^u^ + lusa:) + 0{h^) = 0. (219) 

The term in h disappears when 5 — —2a. Substituting a — — |, ^ = 1, into (217) and (219), 
we get 

ii — (l + ah^u + Ph'^u'^) {|''^-2 — U-i + ui — \u2 + \ah?' \u-1U-2 + + u{u-\ — u\) 

-uj - 1*11*2] + [u^_i{u-2 + u)- uj{u + U2)] } , (220) 

and 

Ut + h?{Qauu^ + 6/?m\^ + u^j) = 0, (221) 

where the 0{h'^) terms were ignored. Using a scale, t — > h'^t, we get (189) and (190). 

Remarks. Step 2 is well suited for a CAS. Solving the system, as outlined in Step 3, is a 
challenging task, in particular, if attempted with pen and paper. The system consists of two 
DDEs and eighteen difference equations. None of the CAS has a build-in solver for such mixed 
systems. A fully automated solution is therefore impossible. We used a feedback mechanism 
which mimics what one would do by hand: Solve the simplest difference equations; enter that 
partial solution; let CAS simplify the entire system; repeat the process until all difference 
equations are satisfied and the two DDEs are simplified as far as possible. Continuing with 
(214) and (215), steps 4 and 5 were straightforward to implement. The five steps have been 
implemented in Mathematica [50]. Starting from (191), the code generates the Gardner lattice 
(220). The code could be modified to assist in the derivation of other completely integrable 
DDEs, such as discrete versions of the nonlinear Schrodinger and sine-Gordon equations. 
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13.2 Dilation Invariance of the Gardner Lattice 

Since (189) is not uniform in rank we must introduce auxiliary parameters with weight. This 
can be done in several ways. One of the possibihties is to replace (189) by 

u = (j + au + (3u^) ^ ^^'"-2 ~ + ^1 ~ ^""2^ + ^ [u-lU-2 + u\ + u{u-i - Ui) 



-ul - U1U2] + ^ [u\{u_2 + U)- ul{u + M2)] I , 



(222) 



where we have set h = 1 (by scaling) and introduced a parameter 7. 

Expressing uniformity of rank, setting W{Dt) = 1, and solving the linear system for the 
weights, one finds that W{u) = W{a) = |, W{(3) = 0, and W{'y) = |. So, we do not need a 
scale on (3 and (222) is invariant under the scaling symmetry 

(t, u, a, 7) {X~H, A3m, A^a, A^). (223) 

For /3 = 0, (222) reduces to 

It = (7 + au) I7 Qm-2 - + ui- ^U2^ + ^ [it-iit-2 + u^-i + u{u^i - Ui) 

-ul - U1U2] } , (224) 

which is a completely integrable discretization of the KdV equation, Ut + GauUx + u-^^ = 0. 
Computing the weights, one can set W{a) = 0, which leads to W{u) = W{-f) = \. So, (224) 

is invariant under the scaling symmetry (t, m, 7) — >• (A~^t, A^m, A^7). For a = 0, (222) reduces 
to 

it = (7 + (5u^) I7 ^""-2 - u-i + - ^^2^ + ^ [u\{u-2 + m) - ul{u + M2)] I , (225) 

which is a completely integrable discretization of the mKdV equation, Ut + Q(5v?Ux + u^x — 0. 
In this case, one can set W{(3) — 0. Then W{u) = ^ and W{'y) — |. Thus, (225) is invariant 
under the scahng symmetry (i, u, 7) — > {X~^t, X^u, A27). 



13.3 Conservation Laws of the Gardner Lattice 

One can either apply the method of Section 12 directly to (189) or, alternatively, apply to 
(222) the technique based on dilation invariance outlined in Sections 9.2, 11.2, and 11.4. In 
particular, one can use the "divide and conquer" strategy of Section 11.5 to split candidate 
densities into smaller pieces. Computational details can be found in [34] The results below 
were computed [58] with the method in Section 12. For q — shifts there are two (non- 
polynomial) density-flux pairs: 

pf ^ = In (1 + ah^ u + f3h^ u^), (226) 

jf'^ = ~| 1^ ('"-2 - ""-i - U + Ui) + J- {2U-2U - 4U-iU + 2U-iUi) 

-\-a^h {u-i{u-2 + U-i -\-u) -\- u{u-i + u + ui)) + aP h (u\u-2 + 2u-2U-iu + Su'^^u 
+3u-iU^ + 2u-iUUi + u^Uij + 2P^ hu-iu(u^2U'-i + U-iU + uuij (227) 
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and 

(228) 



pl^^ = arctanh ( M*^ + 2/3tt)^ 

^2°^ = — 4/^1^ (2m_i — M_2 — 2mi + U2) + a (m_2M-i + u\ + 2m_im + + uuij 

+/3 (ii^i(ii-2 + li) + 1*^(1^-1 + ui)) }. (229) 
The next two (of infinitely many polynomial) densities are 

^ uui + ^ u, (230) 

p(^) = UU2 (1 + ah^ui + f^h'^ul) + ah?uui{u + -ui) + \ph^v?u\ 

where the associated fluxes have been omitted due to length. 
Special Cases. 

We consider two important special cases. The first few densities for (189) with /3 = are: 

pf^ ^ \n{l + ah^u), pf=u, (232) 
p^^^ = ^u'^ + uui, p^"^^ = UU2{1 + ah'^ui) + ah^u{ul + uui + ^u"^). (233) 

The first few densities for (189) with a — are 

pS°^ = ln(l + /?/iV), pf' = arctan(v^/iM), (234) 
p(^) = -UMi, p(2) = iiM2(l + /5/i^M?) + I (^235) 

14 Additional Examples of Nonlinear DDEs 

14.1 The Bogoyavlenskii Lattice 

The Bogoyavlenskii lattice [21] and [90, Eq. (17.1.2)], 

(p p \ 

Y[uj-Y[u-j\, (236) 
j=i j=i / 

is a generalization of the KvM lattice (110). For p — 2, lattice (236) becomes 

ii = u{uiU2 — U-1U-2), (237) 
which is invariant under the following scaling symmetry 

{t,u) ^ {X-H,X^u). (238) 
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Lattice (237) has the following density-flux pairs (of inflnitely many): 

= \nu, J(°) = -(m_iM_2 + + ««!), (239) 

= u, = -'um_i(m_2 + mi), (240) 

p(^) — uui, J^"^^ — —U-iuui{u-2 + u + U2), (241) 

= ^^^^^2 + M2W3), (242) 

+ulu3 + U2U3U4). (243) 



For (237), we also computed the densities p^^^ through p^^\ Every time the rank increases by 
one, the number of terms in the density increases by a factor three. For example, p^^^ has 
2187 terms and the highest shift is 15. 

14.2 The Belov-Chaltikian Lattice 

The Belov-Chahikian lattice [18, Eq. (12)], 

ii — u{ui — U-i) + v^i — V, V — v{u2 — u^i), (244) 

in invariant under the scaling symmetry 

{t, u, v) {X-H, Xu, Xh) . (245) 

The first few density-fiux pairs (of infinitely many) are 

p(^) = u, = -M.iM + f;.!, (246) 

p^^^ = ^u^ + uui — V, J^^-* = —U-iu^ — U-iuui + uv-i + uiV-i + U-iv, (247) 
— u(^u'^ + uui + ul + U1U2 — V-2 — v^i — V — vi), (248) 

where J^^^ has been omitted due to length. Our results match these in [84] . 

14.3 The Blaszak-Marciniak Lattices 

In [19], Blaszak and Marciniak used the R matrix approach to derive families of integrable 
lattices involving three and four fields. Below we consider two cases involving three fields. 
Examples based on four fields could be dealt with in a similar fashion [113]. 
The Blaszak-Marciniak three field lattice I [84, Eq. (2)], 

M = Wi — V = U-iW-i — uw, w = w{v — Vi), (249) 

is invariant under the scaling symmetry 

{t, u, V, w) {X~^t, X^u, Xv, X^w). (250) 
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J« = 


—W-i — w, 






J(2) = 


U-iW-1, 




— ^v^ + uw, 




U-iVW-1 — w 




— ^V^ + UVW + UViW - 


- WWi, 


J{4) 


= W_i('U_it'^ 




— vw — Viw). 



We computed the following density- flux pairs of (249), which is a completely integrable lattice: 

(251) 
(252) 
(253) 
(254) 

(255) 
(256) 

Our results confirm those in [84] and [112]. The Blaszak-Marciniak three field lattice II [112, 
Eq. (1.4)], 

it = vi — V + u{w^i — w), V — v{w-2 — w), w — ui — u, (257) 
is invariant under the scaling symmetry 

{t, u, V, w) {X-H, X^u, AS, Xw). (258) 

The first few density-flux pairs for (257), which is completely integrable, are 

= w, J(^) = -u, (259) 

= l^y2 _ j(2) ^ ^ _ (^260) 

p(^) — ^w'^ + V — uw^i — uw, J^^^ — U^iU + VW^2 + vw^i — UW^i, (261) 

p^^^ — jw^ + \u^ + uui + vw-2 + vw-\ — uwti + vw — uw-iw — uw^, (262) 
J^^^ — —U-2V — uv — u^iVi + vw'^_2 + 2u-iuw-i + VW-2W-1 + vty^i — uw^_-^ + u^iuw. (263) 

Our results agree with those in [84] and [112]. 

15 Software to Compute Conservation Laws for PDEs 
and DDEs 

We first discuss our packages for conservation laws of PDEs and DDEs, followed by a brief 
summary of symbolic codes developed by other researchers. 

15.1 Our Mathematica and Maple Software 

The package TransPDEDensityFlux.m [9], automates the computation of conservation laws of 
nonlinear PDEs in (1 + 1) dimensions. In addition to polynomial PDEs, the software can han- 
dle PDEs with transcendental nonlinearities. The results in Sections 3 and 5 were computed 
with TransPDEDensityFlux.m and cross-checked with the newest version of condens.m, in- 
troduced in [38]. We used TransPDEDensityFlux.m to compute the density-flux pairs for the 
examples in Sections 6 and 7. Details about the algorithm and a discussion of implementation 
issues can be found in [8]. 

The code DDEDensityFlux.m [35] was used to compute the conservation laws in Sections 9 
and 11. The results were cross-checked with the latest version of diffdens.m, featured in 
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[39]. Using multiple scales, the efficiency of DDEDensityFlux.m was drastically improved. 
Nonetheless, the algorithms [34] within DDEDensityFlux.m are impractical for finding den- 
sities and fluxes of high rank. Therefore, we used the new Maple library discrete [56] to 
compute the results in Sections 13 and 14. 

Some of the features of earlier versions of condens .m and dif f dens .m were combined into 
the InvariantsSymmetries .m [39, 41], which allows one to compute generalized symmetries 
as well as conserved densities (but no fluxes). InvariantsSymmetries .m is available from 
MathSource, the Mathematica program bank of Wolfram Research, Inc. 

Our Mathematica packages and notebooks are available at [48] and Hickman's Maple 
code is available at [56]. Our Mathematica codes for the continuous and discrete Euler and 
homotopy operators in one dimension are available at [49]. We are currently designing a 
comprehensive package to compute conservation laws of PDEs in multiple space dimensions 
[45, 51, 83]. 

Our codes have been used in a variety of research projects. For example, condens. m 
[38] was used by Sakovich and Tsuchida [85, 95] to compute conservation laws of nonlinear 
Schrodingcr equations. In [84], Sahadevan and Khousalya use the algorithms of dif f dens. m 
[42] and InvariantsSymmetries .m [39, 41] to compute conserved densities of the Belov- 
Chaltikian and Blaszak-Marciniak lattices. Ergeng and Karasozen [33] used our software in 
the design of Poisson integrators for Volterra lattices. 

15.2 Software Packages of Other Researchers 

Our Mathematica code for conservation laws of PDEs has been "translated" [108, 109, 110] 
into a Maple package, called CONSLAW, which only handles PDEs in (1 + 1) dimensions. Based 
on the concept of dilation invariance and the method of undetermined coefficients, similar 
software was developed by Dcconinck and Nivala [27] and Yang et al. [107]. Our algorithms 
[38, 42] for DDEs have been adapted to fully-discretized equations [36, 37]. 

There are several algorithms (see e.g. [101]) to symbolically compute conservation laws 
of nonhnear PDEs but few have been fully implemented in CAS. Wolf's package ConLaw 
[101, 102, 106] computes first integrals of ODEs and conservation laws of PDEs. ConLaw 
uses the REDUCE package CRACK [103, 104, 105], which contains tools to solve overdetermined 
systems of PDEs. Wolf's application packages heavily rely on the capabilities of CRACK, 
which took years to develop and perfect. Unfortunately, no such package is available in 
Mathematica. 

A common approach is to use the link between conservation laws and symmetries as stated 
in Noether's theorem [15, 66, 81]. Among the newest software based on that approach is the 
Maple code GeM [25] by Cheviakov [23. 24]. which allows one to compute conservation laws 
of systems of ODEs and PDEs based on the knowledge of generalized symmetries. However, 
the computation of such symmetries [47] is as difficult a task as the direct computation of 
conservation laws for it requires solving systems of overdetermined PDEs with, e.g., the Rif 
package [47, 100] . Some methods circumvent the existence of a variational principle (required 
by Noether's theorem) [12, 20, 101, 106] but they still rely on software to solve ODEs or PDEs. 

The package DE_APPLS [16, 77] also offers commands for constructing conservation laws 
from (variational) symmetries by Noether's theorem, but the computation is not fully auto- 
mated. Likewise, the package Noether [43] in Maple allows one to compute conservation laws 
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from infinitesimal symmetry generators corresponding to (simple) Lagrangians. 

Based on the formal symmetry approach, Sokolov and Shabat [89], Mikhailov et al. 
[74, 75], and Adler et al. [10] classified completely integrable PDEs and DDEs in (1 + 1) 
dimensions. Unfortunately, the software used (see [38]) in the classification is obsolete. 

For completeness, we also mention the packages Jets by Marvan [73] and Vessiot by 
Anderson [16, 77]. Both are general purpose suites of Maple packages for computations on 
jet spaces. The commands within Jets and Vessiot use differential forms and advanced 
concepts from differential geometry. By avoiding differential forms, our codes were readily 
adaptable to nonhnear DDEs (not covered in Jets and Vessiot). 

Finally, Deconinck and Nivala [26] developed Maple software for the continuous and dis- 
crete homotopy operators. Their code is available at [28]. 

16 Summary 

We presented methods to symbolically compute conservation laws of nonlinear polynomial 
and transcendental systems of PDEs in (1 + 1) dimensions and polynomial DDEs in one 
discrete variable. 

The first part of this chapter dealt with nonlinear PDEs for which we showed the com- 
putation of densities and fluxes in detail. Using the dilation invariance of the given PDE, 
candidate polynomial densities are constructed as linear combinations with undetermined 
coefficients of scaling invariant building blocks. For polynomial PDEs, the undetermined 
coefficients are constants which must satisfy a linear system of algebraic equations. That 
system will be parameterized by constants appearing in the PDE, if any. For transcendental 
PDEs, the undetermined coefficients are functions which much satisfy a linear system which 
is a mixture of algebraic equations and ODEs. 

The continuous homotopy operator is a powerful, algorithmic tool to compute fluxes ex- 
plicitly. Indeed, the homotopy operator handles integration by parts which allowed us to 
invert the total derivative operator. The methods for polynomial PDEs are illustrated with 
classical examples such as the KdV and Boussinesq equations and the Drinfel'd-Sokolov- 
Wilson system. The computation of conservation laws of system with transcendental nonlin- 
earities is applied to sine-Gordon, sinh-Gordon, and Liouville equations. 

In the second part we dealt with the symbolic computation of conservation laws of non- 
linear DDEs. Again, we used the scaling symmetries of the DDE and the method of un- 
determined coefficients to find densities and fiuxes. In analogy with the continuous case, to 
compute the flux one could use the discrete homotopy operator, which handles summation by 
parts and inverts the forward difference operator. However, in comparison with the "splitting 
and shifting" technique, the discrete Euler and homotopy operators are inefficient tools for the 
symbolic computation of conservation laws of DDEs. The undetermined coefficient method is 
illustrated with classical examples such as the Kac-van Moerbeke, Toda and Ablowitz-Ladik 
lattices. 

There is a fundamental difference between the continuous and discrete cases in the way 
densities (of selected rank) are constructed. The total derivative has a weight but the shift 
operator does not. Consequently, a density of a PDE is bounded in order (with respect to 
x). Unfortunately, there is no a priori bound on the number of shifts in the density, unless 
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a leading order analysis is carried out. To overcome this difficulty and other shortcomings 
of the undetermined coefficient method, we presented a new method to compute conserved 
densities of DDEs. That method no longer uses dilation invariance and is no longer restricted 
to polynomial conservation laws. Instead of building a candidate density with undetermined 
coefficients, one first computes the leading order term in the density and, second, generates 
the correction terms of lower order. The method is fast and efficient since no unnecessary 
terms are computed. The new method was illustrated using a modified Volterra lattice as an 
example, and applied to lattices due to Bogoyavlenskii, Belov-Chaltikian, Blaszak-Marciniak, 
and Gardner. A derivation of the latter lattice was also given. 
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